A Fast and Accurate Failure Frequency Approximation for $k$-Terminal Reliability Systems
Anoosheh Heidarzadeh Alex Sprintson Chanan Singh
cs.DS cs.NI








A Fast and Accurate Failure Frequency Approximation for k-Terminal Reliability Systems
    Anoosheh Heidarzadeh, Member, IEEE,
Alex Sprintson, Senior Member, IEEE, and Chanan Singh Fellow, IEEE
The authors are with the Department of Electrical and Computer Engineering, Texas A&M University, College Station, TX 77843 (E-mail: {anoosheh,spalex,singh}@tamu.edu).
    December 30, 2023
==================================================================================================================================================================================================================================================================================

plain

This paper considers the problem of approximating the failure frequency of large-scale composite k-terminal reliability systems. In such systems, the nodes (k of which are terminals) are connected through components which are subject to random failure and repair processes. At any time, a system failure occurs if the surviving system fails to connect all the k terminals together. We assume that each component's up-times and down-times follow statistically independent stationary random processes, and these processes are statistically independent across the components. In this setting, the exact computation of failure frequency is known to be computationally intractable (NP-hard). In this work, we present an algorithm to approximate the failure frequency for any given multiplicative error factor that runs in polynomial time in the number of (minimal) cutsets. Moreover, for the special case of all-terminal reliability systems, i.e., where all nodes are terminals, we propose an algorithm for approximating the failure frequency within an arbitrary multiplicative error that runs in polynomial time in the number of nodes (which can be much smaller than the number of cutsets). In addition, our simulation results confirm that the proposed method is much faster and more accurate than the Monte Carlo simulation technique for approximating the failure frequency. 

 







k-terminal reliability systems, failure probability and failure frequency, polynomial-time approximation algorithms, minimum cutsets and near-minimum cutsets. 




§ ACRONYMS





DNF     disjunctive normal form

KLM     Karp-Luby-Madras

MCS     Monte Carlo simulation

NP-hard     non-deterministic polynomial-time hard

RGC     recursive generalized contraction




§ NOTATION





Ψ     a composite system

m     number of components

n     number of nodes

k     number of terminal nodes

[i]     set of integers {1,…,i}

λ_i     failure rate of component i

μ_i     repair rate of component i

λ_max     maximum failure rate of a component

μ_min     minimum repair rate of a component

λ     sum of failure rates of all components

μ     sum of repair rates of all components

p_i     probability of component i being unavailable

w_i     weight of component i

w_max     maximum weight of a component

w     sum of weights of all components





















C     set of all minimal cutsets

N     number of minimal cutsets

𝒞_j     jth minimal cutset

s^*     minimum size of a cutset

w(𝒞)     weight of a cutset 𝒞

w^*     minimum weight of a cutset

p^*     maximum failure probability of a cutset

p(𝒞_ℐ)     probability of all components in the collection of cutsets {𝒞_j}_j∈ℐ being unavailable


ϵ     target approximation error factor


δ     target approximation error probability

P_f     probability of system failure in steady-state

P_f^+, P_f^-     first-order upper- and lower-bound on P_f

P̂_f     an approximation of P_f using bounding technique



P̃_f     an approximation of P_f

F_f     frequency of system failure in steady-state

F_f^+, F_f^-     first-order upper- and lower-bound on F_f

F̂_f     an approximation of F_f using bounding technique



F̃_f     an approximation of F_f

P     probability of all components of some cutset being unavailable and unexposed

P̃     an approximation of P

C^(α)     set of all α-min cutsets for arbitrary α≥ 1

N^(α)     number of α-min cutsets

P^(α)     probability of all components of some α-min cutset being unavailable and unexposed 

P_f^(α)     probability of all components of some α-min cutset being unavailable

P̃_f^MC     an approximation of P_f using Monte Carlo simulation

F̃_f^MC     an approximation of F_f using Monte Carlo simulation

s     (Boolean) system-state vector

𝒮_f     set of all system-states in which the system is unavailable

ℐ(s)     set of unavailable components in system-state s

p(s)     probability of system-state s

poly(N)     a polynomial function in N

poly(n)     a polynomial function in n





§ NOMENCLATURE





α-min cutset     a minimal cutset of weight no greater than α times the minimum cutset weight

(ϵ,δ)-approx.     a multiplicative approximation with error factor at most ϵ and error probability at most δ
	













§ INTRODUCTION

Consider a composite system whose n nodes, consisting of k (2≤ k≤ n) terminals and n-k relays, are connected through components that are subject to statistically-independent continuous-time stationary failure/repair random processes. At any given time, each component is either operational or not, and the system fails if the surviving system of operational components does not connect all terminals. The probability of failure (P_f) and the frequency of failure (F_f) are two important measures of reliability of such systems <cit.>. These quantities are very useful to derive other reliability measures such as mean down-time and mean cycle-time <cit.>. 




Numerous algorithms were previously designed for computing P_f <cit.> and F_f <cit.>. The proposed algorithms, however, become intractable in large-scale systems because the computational complexity grows very quickly with the number of nodes in the system (i.e., n). Specifically, the exact computation of these quantities was shown to be NP-hard <cit.>. To overcome this challenge, various methods were developed to approximate P_f and F_f. 



The existing techniques for approximating P_f are based on the Monte Carlo simulation (MCS) <cit.> and rare-events simulation <cit.>. Similar methods were used to approximate F_f in <cit.> and <cit.>. Notwithstanding, the only existing technique which can provably approximate P_f within an arbitrary multiplicative error and runs in polynomial time in the number of cutsets (i.e., minimal collections of components whose failure results in a loss of connectivity of some terminals from the rest) was proposed in <cit.>. To our knowledge, no such computationally efficient algorithm with provable guarantees was previously proposed for approximating F_f with an arbitrary multiplicative error factor. 





The algorithm of <cit.> estimates P_f for the settings where the failure/repair random process for each component is stationary. In such settings, each component is available or unavailable at any time instant, independently from other components, with some constant probability (independent of time). Then, P_f is equal to the probability that all components in a cutset are unavailable. The algorithm of <cit.> relies on the fact that P_f can be written as the truth probability of some disjunctive normal form (DNF) formula, and a multiplicative approximation of this probability can be computed in polynomial time in the number of cutsets in the system. However, the number of cutsets can be exponential in n, and the application of this technique is not practical in such cases. 







Interestingly, for the special case where all n nodes are terminals (i.e., k=n), it was shown <cit.> that a multiplicative approximation of P_f can be computed in polynomial time in n. (Such results do not exist for more general cases where only a subset of nodes are terminals.) One example of such all-terminal reliability systems is the electricity distribution networks where the terminals represent either a supply point or a major load point. The feeders between these terminals are the components that can fail and be repaired. The loss of supply to any major load point is considered a system failure. At the distribution level there may be only one feed point for many of such networks, and in such cases, the networks are all-terminal reliability systems. (For a network with more than one feed point, depending on its topology, the network may or may not be an all-terminal reliability system.) Other examples of such systems can be found in communication networks, computer networks, and transportation networks <cit.>.

The main ideas in <cit.> can be summarized as follows: (i) the number of weak cutsets in a system, i.e., those cutsets with higher probability of failure, is polynomial in n, and the enumeration of all such cutsets can be done in polynomial time in n, and (ii) the probability that all components of a weak cutset are unavailable provides a multiplicative approximation of P_f. The idea of using weak cutsets was also used recently in <cit.> for computing a bounding-type approximation of P_f. Nevertheless, there is no apparent connection between F_f and the truth probability of a DNF formula. This implies the need for a novel technique for approximating F_f.

The main contributions of this work are as follows:


  * We present a new algorithm that runs in polynomial time in the number of cutsets in the system, and approximates F_f for k-terminal reliability systems (for any k) within any given multiplicative error. Moreover, for the special case where all nodes are terminals, we propose a new algorithm to compute an approximation of F_f for any given multiplicative error factor that runs in polynomial time in n. To our knowledge, this is the first and only fast (polynomial time in n) and accurate (with arbitrary provable guarantees) algorithm for approximating F_f for all-terminal reliability systems. 

  * We adapt the machinery of <cit.>, which was tailored to approximate P_f, in a non-trivial way to approximate F_f. In particular, we present a new transformation to obtain F_f from P_f by introducing an auxiliary probability P, and re-writing F_f as a scalar multiple of the difference between P_f and P. To the best of our knowledge, this connection was not previously reported in the literature. This transformation enables us to convert the problem of approximating F_f to the two sub-problems of approximating P_f and P, each of which can be linked to a DNF formula via a carefully designed random process. 

  * It is well known that P_f can be thought of as the probability that, under a random sampling process, all components of a cutset in the system are unavailable. To relate P to a random process, we define a new random process, referred to as the exposure process, statistically independent from the sampling process, such that P can be thought of as the probability that, under the sampling and exposure processes simultaneously, all components of a cutset in the system are unavailable and unexposed. 

  * We identify the weak cutsets of a system that satisfy the following property: under the sampling and exposure processes simultaneously, the probability that all components of a weak cutset are unavailable and unexposed provides multiplicative approximations of P_f and P, and subsequently, a multiplicative approximation of F_f for any arbitrary error factor. We also prove that the set of weak cutsets required for approximating F_f contains (and is greater, but not more than a factor polynomial in n, than) the set of weak cutsets previously identified in <cit.> for approximating P_f.

  * We show that MCS can provide an additive approximation of F_f in polynomial time; whereas, using MCS, F_f cannot be approximated in polynomial time within a multiplicative error. This suggests that approximating F_f with a multiplicative error factor, when compared to approximating F_f with an additive error factor, is computationally more expensive, and hence more challenging.

  * We compare the proposed technique and the MCS technique via simulations for approximating F_f  for a 3× 3 (9 nodes) grid network and for the layer 3 (20 nodes) and the layer 2 (35 nodes) of the Internet2 network <cit.>. (The results of the bounding technique are also given for reference.) Our simulation results show that the proposed technique, when compared to the MCS technique, provides approximations with higher accuracy, for less running time. 


The rest of the paper is organized as follows. Section <ref> gives basic definitions/notations and the problem formulation. In Section <ref>, we overview the concepts of failure and repair rates, and cutsets and near-minimum cutsets. Sections <ref> and <ref> discuss the previous works and the proposed algorithms for k-terminal and all-terminal reliability systems, respectively. In Section <ref>, we present our simulations results and compare the proposed technique for all-terminal reliability systems with the Monte Carlo simulation technique. Section <ref> concludes the paper and discusses some open problems. The proofs of some lemmas are deferred to the appendix. 
	


§ PROBLEM SETUP

Let Ψ be a system with n nodes, k (2≤ k≤ n) of which are terminals and the rest are relays, and m components (i.e., the edges connecting the nodes). The system Ψ is called a k-terminal reliability system, and for the case of k=n, the system Ψ is called an all-terminal reliability system. Fig. <ref> depicts a 3× 3 grid network as an example of an all-terminal reliability system with n=k=9 and m=12. (Note that the regularity of the system in Fig. <ref> is not a requirement for the applicability of the proposed techniques in this paper.)

Let [m]≜{1,…,m} be the index set of components in the system Ψ. Each component i∈[m] is assumed to have two states: available (up) and unavaiable (down). We consider the setting that each component of the system Ψ is subject to a statistically-independent continuous-time failure/repair random process described shortly, and study the steady-state behavior of the system Ψ under such random processes.





 §.§ Failure/Repair Random Processes

We assume that every component is initially available. Note that for a stationary ergodic stochastic process, such as a two-state process, the (steady-state) failure probability and (steady-state) failure frequency do not depend on the initial conditions <cit.>. As time evolves, each component i becomes unavailable after a random period of time, distributed arbitrarily with mean 1/λ_i (for arbitrary λ_i>0), and it becomes available again after a random period of time, distributed arbitrarily with mean 1/μ_i (for arbitrary μ_i>0). This process, for every component, continues over time, statistically independent from other components. 

Note that we do not restrict the distributions of up-times and down-times to be exponential. As shown in <cit.>, the (steady-state) failure probability and (steady-state) failure frequency of a system with two-state components depend only on the mean-up-times {1/λ_i} and mean-down-times {1/μ_i} of the components in the system, and not on the up/down-time distribution functions per se. 








 §.§ Failure Probability and Failure Frequency

At any time, the (sub-) system of Ψ including all n nodes restricted only to available components (i.e., the original system excluding unavailable components), is referred to as the surviving system. At any time, the system is unavailable if the surviving system fails to connect all k terminals. The (steady-state) probability of failure is the probability that the system is unavailable, and the (steady-state) frequency of failure is the expected number of times per unit time (i.e., the expected rate) that the system becomes unavailable <cit.>.



 §.§ Problem Statement

For arbitrary ϵ>0 and 0<δ<1, we consider the problem of computing (ϵ,δ)-approximations of P_f and F_f (with an error factor of at most ϵ and an error probability of at most δ), denoted by P̃_f and F̃_f, respectively, defined as 
    {|P̃_f-P_f|≥ϵ P_f }≤δ,
 and 
    {|F̃_f-F_f|≥ϵ F_f }≤δ.
 



 §.§ Complexity Notation

Throughout the paper, we follow the conventional use of notations in complexity theory as described below. 

For two arbitrary functions f(x) and g(x) of variable x, we write: (i) f(x)=O(g(x)) if for sufficiently large x, |f(x)|≤ c· |g(x)| for some constant c>0; (ii) f(x)=Ω(g(x)) if for sufficiently large x, f(x)≥ c· g(x) for some constant c>0; (iii) f(x)=Θ(g(x)) if f(x)=O(g(x)) and f(x)=Ω(g(x)); and (iv) f(x)=o(g(x)) if for any constant c>0, |f(x)|≤ c· |g(x)| for sufficiently large x.



















 


§ PRELIMINARIES



 §.§ Failure/Repair Rates and Unavailability Probabilities


We refer to λ_i and μ_i as the failure rate and the repair rate of component i, respectively. We assume that λ_i and μ_i are independent of time (i.e., the failure/repair random process is stationary). We also assume that λ_i and μ_i do not depend on n. Since λ_i and μ_i are independent of time, the probability that each component i is unavailable at any given time is equal to p_i≜λ_i/(λ_i+μ_i). Note that the probabilities {p_i} are independent of time. We refer to p_i as the unavailability probability of component i, and refer to w_i=-log p_i as the weight of component i, where the symbol “log”, here and throughout the paper, refers to the natural logarithm. Since 0<p_i<1 for all i, the weights {w_i} are all non-negative. 



Define λ_max≜max_i∈ [m]λ_i as the maximum failure rate of a component. Similarly, define μ_min as the minimum repair rate of a component. For simplifying the arguments, we assume that μ_min/λ_max> m-1. Define w_max as the maximum weight of a component. Moreover, let λ≜∑_i=1^mλ_i, μ≜∑_i=1^mμ_i, and w≜∑_i=1^m w_i. 

We assume, without loss of generality, that there are no parallel components in the system. Otherwise, if there exist l parallel components {i_j}_1≤ j≤ l with failure rates {λ_i_j}, repair rates {μ_i_j}, and unavailability probabilities {p_i_j}, we can replace them all by one component with unavailability probability p=∏_j p_i_j, repair rate ∑_jμ_i_j, and failure rate (p∑_jμ_i_j)/(1-p). This transformation does not change the failure probability and the failure frequency of a system <cit.>. Note that this assumption is made only for the ease of exposition, and is not a requirement for the proposed algorithms. 










 §.§ Cutsets and Near-Minimum Cutsets

A set of components in a system is a cutset if the unavailability of those components yields the unavailability of the system. Moreover, a cutset is a minimal cutset if it does not contain any other cutsets. Hereafter, we use the term “cutset” as a shorthand for “minimal cutset.” Let N be the number of cutsets in a system, and let C ≜{𝒞_1,…,𝒞_N} be the set of all cutsets in the system. Let 𝒞_j ≜{i_1,…,i_|𝒞_j|}, where i_1,…,i_|𝒞_j| are the indices of components in 𝒞_j, and |𝒞_j| is the size of 𝒞_j, i.e., the number of components in 𝒞_j. Hereafter, we assume that N>1 (otherwise, computing P_f and F_f is trivial). For example, in Fig. <ref>, there exist 4 and 16 cutsets of sizes 2 and 3, respectively, as enumerated in Table <ref>. 

For any ℐ⊆ [N], let 𝒞_ℐ≜∪_j∈ℐ𝒞_j. The probability of failure of cutsets {𝒞_j}_j∈ℐ, denoted by p(𝒞_ℐ), is equal to the probability that all components in 𝒞_j, for all j∈ℐ, are unavailable, i.e., p(𝒞_ℐ) ≜∏_i∈𝒞_ℐ p_i. Note that p(𝒞_j) is the probability of unavailability of all components in 𝒞_j. Let w(𝒞_j), the weight of 𝒞_j, be the sum of the weights of all components in 𝒞_j, i.e., w(𝒞_j) ≜∑_i∈𝒞_j w_i. Note that w(𝒞_j) = -log(p(𝒞_j)). Let w^*≜min_𝒞∈ C w(𝒞) and p^*≜max_𝒞∈ C p(𝒞) be the minimum weight and maximum failure probability of a cutset, respectively. Let s^*≜min_𝒞∈  C |𝒞| be the minimum size of a cutset. Note that s^*≥min{max{w^*/w_max,1},m-1}. 






For any constant α≥ 1 (independent of n and m), let C^(α)≜{𝒞∈ C: w(𝒞)≤α w^*} be the set of all (minimal) cutsets in a system of weight less than or equal to α w^*. Let N^(α) be the number of cutsets in C^(α). We refer to the cutsets in C^(α) as α-min cutsets. For simplicity, we refer to 1-min cutsets as min-cutsets. Note that w(𝒞)=w^* and p(𝒞)=p^* for all min-cutsets 𝒞∈ C.

For the cases in which all components have the same weight w_0, there is a one-to-one map between the weights and the sizes of the cutsets (i.e., the weight of each cutset is equal to the size of that cutset times w_0). For example, for the case that all components in Fig. <ref> have unit weight (i.e., all cutsets of size 2 or size 3 have weight 2 or weight 3, respectively), Table <ref> enumerates the 1.5-min cutsets of the system in Fig. <ref>. 


















§ A POLY(N)-TIME APPROXIMATION ALGORITHM FOR K-TERMINAL RELIABILITY SYSTEMS




























 §.§ Background

In this section, we overview the inclusion-exclusion based formulas for P_f and F_f, and the bounding technique for approximating P_f and F_f. The reader familiar with these concepts can skip this section. 

By the cutset approach <cit.> based on the inclusion-exclusion principle, P_f and F_f can be written as:

    P_f = ∑_l=1^N((-1)^l+1∑_ℐ⊆ [N]: |ℐ|=l p(𝒞_ℐ))
 and 

    F_f = ∑_l=1^N((-1)^l+1∑_ℐ⊆ [N]: |ℐ|=l(p(𝒞_ℐ)∑_i∈𝒞_ℐμ_i))
 (The details of derivation of the formulas (<ref>) and (<ref>) can be found in <cit.>.) Note that, depending on the topology (e.g., ring topology), the number of cutsets (N) in some systems is only poly(n), and all cutsets can be enumerated in O(N) time <cit.>. However, the complexity of computing P_f and F_f using the formulas (<ref>) and (<ref>) may be still unaffordable. (The number of terms in formulas (<ref>) and (<ref>) is exponential in the number of cutsets (N) and double-exponential in the number of nodes (n). Thus, there does not exist a poly(n)-time (i.e., with running time polynomial in n) nor a poly(N)-time (i.e., with running time polynomial in N) algorithm for computing P_f and F_f from (<ref>) and (<ref>), directly.) For such systems, one may use the bounding technique to approximate P_f or F_f. 

The bounding technique is one of the most common approaches to provide upper and lower bounds on P_f or F_f via truncating the formula (<ref>) or (<ref>), respectively <cit.>. 



Let 

    P_f^+ = ∑_1≤ j≤ N p(𝒞_j),
 and 

    P_f^- = ∑_1≤ j≤ N p(𝒞_j)-∑_1≤ j_1<j_2≤ N p(𝒞_j_1∪𝒞_j_2).
 Similarly, let 

    F_f^+ =∑_1≤ j≤ N(p(𝒞_j) ∑_i∈𝒞_jμ_i),
 and 

    F_f^- = ∑_1≤ j≤ N(p(𝒞_j) ∑_i∈𝒞_jμ_i) - ∑_1≤ j_1<j_2≤ N(p(𝒞_j_1∪𝒞_j_2)∑_i∈𝒞_j_1∪𝒞_j_2μ_i ).

 Then, P_f^+ and P_f^- (or F_f^+ and F_f^-) are the first-order upper-bound and lower-bound on P_f (or F_f), respectively. Using a similar technique by incorporating larger collections of cutsets (instead of singletons or pairs only), one can compute higher-order upper- and lower-bounds on P_f and F_f <cit.>. Such higher-order bounds, when compared to the first-order bounds, are more accurate, but more computationally expensive.

Now a question is what type of guarantee the bounding technique provides on the accuracy of the approximation. To answer this question, let d_P (or d_F) be the maximum number of decimal places up to which P_f^+ and P_f^- (or F_f^+ and F_f^-) match. Let P̂_f = trunc(P_f^+,d_P) = trunc(P_f^-,d_P) and F̂_f=trunc(F_f^+,d_F)=trunc(F_f^-,d_F), where trunc(x,d)=⌊ 10^d· x⌋/10^d, for any real number x≥ 0 and integer d≥ 1. Note that P̂_f and F̂_f are the most accurate estimators of P_f and F_f based on the (first-order) upper and lower bounds. Both bounds P̂_f and F̂_f are computable in poly(N) time, yet neither guarantees an arbitrary approximation error factor. To be more specific, P̂_f and F̂_f are always exact up to d_P and d_F decimal places, respectively, but the approximation error depends on d_P and d_F, and cannot be made arbitrarily small.





 §.§ Main Ideas of the Proposed Algorithm

In this section, we give an overview of the KLM estimator as part of the proposed algorithms, and explain our main ideas. 

Let the symbols “∧” and “∨” denote the logical conjunction (AND) and the logical disjunction (OR), respectively. Let 
    Φ_M=Z_1∨ Z_2∨…∨ Z_M
 be a formula on M>1 variables {Z_j}_j=1^M, where the clause Z_j is a conjunction of literals z_i for some i∈ [m], i.e., Z_j=∧_i∈ I_j z_i for some I_j⊆ [m]. Each literal z_i is either a Boolean variable or the negation of a Boolean variable, and it takes two values: “true” and “false.” The formula Φ_M of such form is said to have disjunctive normal form (DNF). Let P_Z(j) be the probability that the clause Z_j is true (i.e., the literals {z_i}_i∈ I_j are all true), and let P_Z denote the vector [P_Z(1),…,P_Z(M)]. Note that the truth probability of Φ_M, denoted by Π_M, cannot be computed in poly(M)-time (i.e., polynomial-time in M) <cit.>. However, using an unbiased estimator, referred to as KLM, due to Karp, Luby, and Madras <cit.>, one can compute a (ξ,δ)-approximation of Π_M in poly(M) time (for any given ξ>0 and 0<δ<1). In the following, we describe a simple, yet powerful, extension of the KLM estimator. 

The KLM estimator with inputs (Φ_M,P_Z; ξ,δ) proceeds in steps as follows: 


  0. Initialize the counters s and t by setting s=t=1;

  1. Choose a random clause Z_j, with probability of selecting Z_j being equal to P_Z(j)/Q_Z, where Q_Z=∑_j P_Z(j);

  2. Choose a random assignment z={z_1,…,z_m} satisfying the clause Z_j (i.e., z_i is true for all i∈ I_j); 

  3. Compute π_s,t=Q_Z/N(z), where N(z) is the number of clauses that the assignment z satisfies; 

  4. s← s+1  

  5. Repeat Steps 1-4 S=⌈4(M-1)/ξ^2⌉ times;

  6. Compute the mean π_t=(∑_s=1^Sπ_s,t)/S; 

  7. t← t+1 

  8. Repeat Steps 1-7 T = ⌈ 12 log(1/δ)⌉ times;

  9. Return the median of {π_t}_t=1^T.
 

The running time of an obvious implementation of the KLM estimator is O((M^2 m/ξ^2) log(1/δ)). In particular, Step 1 takes O(M) time to run. Step 2 can be simply run in O(m) time; Step 3 takes O(M m) time to run; and each of Steps 1-3 is run ST = O((M/ξ^2)log(1/δ)) times; Step 6 can be run in O(S) time, and this step is run T times, and Step 9 can be run in O(T) time, and this step is run only once. Note that a more sophisticated implementation of the KLM estimator, referred to as self-adjusting, can be run in O((M m/ξ^2)log(1/δ)) time (see, for more details, <cit.>).


Let Π̃_M be the output of the KLM estimator. Then, the following result holds.


Π̃_M is a (ξ,δ)-approximation of Π_M.


The proof follows from similar arguments as those in <cit.>, and can be found in the appendix.	


As can be seen in (<ref>), P_f is the probability of union of a set of events, and consequently, it can be thought of as the probability of satisfying a DNF formula with random Boolean variables. Thus, the KLM estimator can compute an (ϵ,δ)-approximation of P_f for k-terminal reliability systems (for any k) in poly(N) time. However, as one can see in (<ref>), F_f is not the probability of union of any set of events. Thus, F_f cannot be thought of as the probability of satisfying a DNF formula, and the KLM estimator is not directly applicable. This poses a challenge for approximating F_f. To overcome this challenge, we define an auxiliary term P that satisfies the following requirements: (i) P can be interpreted as the probability of the union of an auxiliary set of events, and (ii) F_f is a scalar multiple of the difference between P_f and P. Note that, for the first time in the literature, this work presents a linear connection between P_f and F_f. 

Each of P_f and P is the probability of union of a set of events, and thus can be approximated by the KLM estimator within an arbitrary multiplicative error. Now the main questions are whether P_f and P can be approximated in polynomial time in the number of cutsets, and whether combining such approximations of P_f and P, a multiplicative approximation of F_f with an arbitrary error factor can be computed. In this section, we answer these questions in the affirmative, and propose a poly(N)-time algorithm that computes an (ϵ,δ)-approximation of F_f for k-terminal reliability systems (for any k).  

We define P that satisfies the requirements (i) and (ii) as follows:


    P = ∑_l=1^N((-1)^l+1∑_ℐ⊆ [N]: |ℐ|=l(P(𝒞_ℐ) (1-∑_i∈𝒞_ℐμ_i/μ))).
 By using (<ref>)–(<ref>), it is easy to verify that 

    F_f=(P_f-P)μ,

where F_f and P_f are given by (<ref>) and (<ref>), respectively. Thus, F_f is a scalar multiple of the difference between P_f and P, as desired. Note that P_f> P since F_f>0. (We notice, without proof, that F_f=0 only in a system with p_i=0 for all i, and this contradicts with the assumption that 0<p_i<1.)

Now, we define an auxiliary set of events such that P, defined in (<ref>), is the probability of union of these events. In the steady state, at any given time, each component i is unavailable or available, with probability p_i or 1-p_i, respectively, independent of time <cit.>. Thus, the random process under consideration (Section <ref>) is equivalent to the following one-shot random process over the system Ψ. Each component i, statistically independently from other components, is set to be “unavailable” with probability p_i, and it is set to be “available” otherwise. We refer to this process as sampling. 

We further introduce an auxiliary one-shot random process over the system Ψ as follows. Each component is assumed to have two states: exposed and unexposed. One component, say i, is chosen with probability μ_i/μ, and is set to be exposed; and the rest of the components are set to be unexposed. We refer to this process as exposure. Note that the sampling and exposure processes are statistically independent. 

The intuition behind the sampling and exposure processes is as follows. Each term in P corresponds to a collection of cutsets and expresses the probability that all components in this collection of cutsets are unavailable and unexposed. Moreover, as can be seen in (<ref>), P has the structure of an inclusion-exclusion formula. Thus, it should not be hard to see that P is equal to the probability that all components in some cutset of the system Ψ (under the sampling and exposure processes) are unavailable and unexposed (due to the statistical independence of these processes). 





 §.§ Proposed Algorithm

For a k-terminal reliability system (for arbitrary k), the inputs of the proposed algorithm are the failure rates {λ_i}, the repair rates {μ_i}, and the approximation parameters ϵ>0 and 0<δ<1. The algorithm proceeds in steps as follows:


  0. Initialization:


  0.1 Enumerate all N (minimal) cutsets of the system;

  0.2 Compute s^*, μ_min, λ_max, μ, and ρ=μ_mins^*-λ_max (m-s^*);

  0.3 Take ξ=(ϵ/2)(ρ/μ)


  1. Compute a (ξ,δ/2)-approximation P̃ of P using the KLM estimator; 

  2. Compute a (ξ,δ/2)-approximation P̃_f of P_f using the KLM estimator; 

  3. Return F̃_f = (P̃_f-P̃)μ. 	




The running time of this algorithm is O((N m^3log N) (1/ϵ^2)log(1/δ)). (The running time of each step of the algorithm is given in Section <ref>.)

The details of the computations of P̃_f and P̃ (Steps 1 and 2) are as follows. 



  §.§.§ Computation of P̃

Construct a DNF formula as 
    Φ_N≜ Z_1∨ Z_2∨…∨ Z_N,
 where the clause Z_j, j∈ [N], is the conjunction of two sub-clauses X_j and Y_j (i.e., Z_j=X_j ∧ Y_j). Define 
    X_j≜ x_i_1∧ x_i_2∧…∧ x_i_|𝒞_j|,
 and 
    Y_j≜ y_i_1∧ y_i_2∧…∧ y_i_|𝒞_j|,
 where i_1,…,i_|𝒞_j| are the labels of the components in the cutset 𝒞_j, and the literals x and y are Boolean (random) variables defined as follows. For any (random) assignment of literals, x_i is true with probability p_i=λ_i/(λ_i+μ_i), and it is false, otherwise; and y_i is true for all i, except for one and only one i being chosen with probability μ_i/μ. Define 
    P_Z(j) ≜ p(𝒞_j) (1-∑_i∈𝒞_jμ_i/μ).
 Run the KLM estimator with inputs (Φ_N,P_Z;ξ,δ/2), and denote by P̃ the output. 




  §.§.§ Computation of P̃_f

The method of computing P̃_f is similar to that of P̃, except that a slightly different DNF formula is required. This technique was previously used in <cit.>. Construct a DNF formula as 
    Φ_N,f≜ X_1∨ X_2∨…∨ X_N,
 where the clause X_j, j∈ [N], is defined as before. Define 
    P_X(j)≜ p(𝒞_j).
 Run the KLM estimator with inputs (Φ_N,f,P_X;ξ,δ/2), and denote by P̃_f the output. 



 §.§ Theoretical Analysis



The output of the algorithm in Section <ref>, F̃_f, is an (ϵ,δ)-approximation of F_f. 


 By the definitions of X_j and Y_j as above, it is not hard to see that (i) P_X(j) is the probability that a random assignment x={x_1,…,x_m} satisfies X_j (i.e., X_j is true), and (ii) P_Y(j)≜ 1-∑_i∈𝒞_jμ_i/μ, is the probability that a random assignment y={y_1,…,y_m} satisfies Y_j (i.e., Y_j is true). By (i) and (ii), it follows that a random assignment z=(x,y) satisfies Z_j (i.e., Z_j is true) with probability P_Z(j), defined earlier, since X_j and Y_j are statistically independent, i.e., P_Z(j) = P_X(j)· P_Y(j).  

The formula Φ_N is true so long as some clause Z is true, and thus by the inclusion-exclusion principle, it immediately follows that the truth probability Π_N of Φ_N is equal to P. Thus, in order to approximate P, it suffices to approximate Π_N. Since Φ_N is a DNF formula, Π_N can be approximated using the KLM estimator. By Lemma <ref>, the output of the KLM estimator with inputs (Φ_N,P_Z;ξ,δ/2), denoted by P̃, is a (ξ,δ/2)-approximation of Π_N, or equivalently, P. This yields the following result. 


P̃ is a (ξ,δ/2)-approximation of P.


Moreover, the DNF formula Φ_N,f is true so long as X_j is true for some j. Thus, the probability that Φ_N,f is true, denoted by Π_N,f, is equal to P_f. By a similar argument as above, the following result is immediate.


P̃_f is a (ξ,δ/2)-approximation of P_f.


The results of Lemma <ref> and Lemma <ref> yield 

    {|P̃-P|≥ξ P}≤δ/2,
 and

    {|P̃_f-P_f|≥ξ P_f }≤δ/2.
 By combining (<ref>) and (<ref>), we get 

    {|(P̃_f-P̃)-(P_f-P)|≥ξ (P_f+P) }≤δ.
 Since F_f=(P_f-P)μ (by (<ref>)) and F̃_f = (P̃_f-P̃)μ (by the algorithm), we get 

    {|F̃_f-F_f|≥ξ (P_f+P)μ}≤δ,
 or equivalently, 

    {|F̃_f-F_f|≥ϵ F_f }≤δ,
 so long as 

    ϵ≥ξ(P_f+P/F_f)μ.
 Thus, it suffices to show that our choice of ξ meets the condition (<ref>). To do so, we need to establish lower- and  upper bounds on F_f. This can be done by using the system-state based formulas for P_f and F_f <cit.>. 

Let s≜{s_1,…,s_m} be the state of a system at a given time, where s_i is “true” if the component i is available at that time, and s_i is “false” otherwise. Let 𝒮_f be the set of all states s in which the system is unavailable. Moreover, let p(s) be the probability of the state s, i.e., 

    p(s) ≜∏_i∈ℐ(s) p_i ·∏_i∈ [m]∖ℐ(s) (1-p_i),
 where ℐ(s) and [m]∖ℐ(s) are the set of unavailable and available components, respectively, in the state s. It was shown in <cit.> that P_f and F_f can be written as: 

    P_f = ∑_s∈𝒮_f p(s),
 and 

    F_f = ∑_s∈𝒮_f p(s)( ∑_i∈ℐ(s)μ_i -∑_i∈ [m]∖ℐ(s)λ_i ).
 





Note that |ℐ(s)|≥ s^* for any s∈𝒮_f since (i) all components of at least one cutset are unavailable, and (ii) the size of any cutset is bounded from below by s^*; and |ℐ(s)|≤ m for all s. Thus, 

    μ_mins^*≤∑_i∈ℐ(s)μ_i≤μ
 and 

    0≤∑_i∈ [m]∖ℐ(s)λ_i≤ (m-s^*)λ_max
 for all s∈𝒮_f. By the choice of the algorithm, ρ = μ_mins^*-λ_max (m-s^*) (Step 0.2). Note that ρ = O(m) and ρ=Ω(1) since 1≤ s^*≤ m, and μ_min=Θ(1) and λ_max=Θ(1). Putting (<ref>) and (<ref>) together, it follows that ∑_i∈ℐ(s)μ_i-∑_i∈ [m]∖ℐ(s)λ_i, for all s∈𝒮_f, is lower and upper bounded by ρ and μ, respectively. By (<ref>) and (<ref>), it is then easy to see that 
    P_fρ≤ F_f≤ P_f μ.
 












Using the bounds on F_f in (<ref>), we proceed with the rest of the proof as follows. 

Since F_f≥ P_fρ (by (<ref>)), and P_f+P< 2P_f (by the fact that F_f>0, and so, P_f>P (by (<ref>))), one can see that

    ξ(P_f+P/F_f)μ< 2ξ(μ/ρ).
 Taking ξ=(ϵ/2) (ρ/μ) as in the algorithm, it is obvious that (<ref>), and consequently, (<ref>) hold (by (<ref>)). Thus, F̃_f is an (ϵ,δ)-approximation of F_f.

Note that the smaller is ξ, the larger is the running time of the proposed algorithm. However, the lower bound on F_f and the upper bound on P_f+P dictate the choice of ξ. Thus, the closer are these bounds to the actual values, the more efficient is the computation of the approximation. However, the bounds in (<ref>) follow from a worst-case analysis. For improving on these bounds, the trick is to run the proposed algorithm multiple times, and amplify ξ in each run based on the results of the previous runs. Similar idea was previously used in <cit.>. This is, however, beyond the scope of this paper, and hence not discussed here. 









 §.§ Computational Complexity

The initialization (Step 0) can be run in O(Nlog N) time. In Step 0.1, as was shown in <cit.>, all the N minimal cutsets can be enumerated in O(N) time. In Step 0.2, s^* can be computed in O(Nlog N) time (via sorting), μ_min and λ_max in O(mlog m) time (via sorting), and μ and ρ in O(m) time and O(1) time, respectively.

By the choice of ξ = (ϵ/2)(ρ/μ) in the algorithm, ξ=Ω(ϵ/m) since ρ = Ω(1) and μ = O(m). In Steps 1 and 2, the estimates P̃ and P̃_f can be computed in poly(N) time by running the (self-adjusting) KLM estimator in O((N m/ξ^2)log(1/δ)) time <cit.>, or equivalently, O((N m^3) (1/ϵ^2)log(1/δ)) time. Putting everything together, the proposed algorithm runs in poly(N) time. This technique is useful for systems with poly(n) cutsets. In such systems, F_f can be approximated accurately and efficiently by the proposed algorithm; whereas the accuracy of approximating F_f using the bounding technique may be inadequate (see Section <ref>).





§ A POLY(N)-TIME APPROXIMATION ALGORITHM FOR ALL-TERMINAL RELIABILITY SYSTEMS




 §.§ Background

In this section, we overview the Monte Carlo simulation for approximating P_f and F_f as part of the proposed algorithm. The reader familiar with this concept can skip this section. 




In many systems, there exist exponentially many cutsets and failure states, and it is not practical to identify and enumerate all cutsets or failure states. Thus P_f and F_f cannot be computed from (<ref>) and (<ref>) or (<ref>) and (<ref>), directly. The Monte Carlo simulation (MCS) is useful to approximate P_f and F_f in such cases <cit.>. This technique performs MCS over the system state space (see, for details, (<ref>)–(<ref>)). 

The inputs of the MCS algorithm are the failure rates {λ_i}, the repair rates {μ_i}, and the number of simulation runs, S and T. The MCS algorithm proceeds in steps as follows (see Section <ref> for the definitions of s, p(s), and ℐ(s)): 



  0. Initialize the counters s and t by setting s=t=1;

  1. Choose a random state s, with probability of selecting s being equal to p(s);

  2. If the system is unavailable in the state s, let π_s,t = 1 and φ_s,t=∑_i∈ℐ(s)μ_i-∑_i∈ [m]∖ℐ(s)λ_i; otherwise, let π_s,t =0 and φ_s,t=0; 

  3. s← s+1  

  4. Repeat Steps 1-3 S times;

  5. Compute the means π_t=(∑_s=1^Sπ_s,t)/S and φ_t=(∑_s=1^Sφ_s,t)/S; 

  6. t← t+1

  7. Repeat Steps 1-6 T times;

  8. Return the median of {π_t}_t=1^T and the median of {φ_t}_t=1^T. 
 

The running time of MCS is O(ST(m+n)). In particular, Step 1 can be run in O(m) time; Step 2 can be run in O(m+n) time, e.g., using the breadth first search or the depth first search. Each of Steps 1 and 2 is run ST times; Step 5 can be run in O(S) time, and this step is run T times; Step 6 can be run in O(T) time, and this step is run only once. 

We will show in Section <ref> that MCS can provide approximations of P_f and F_f with arbitrary additive error factors in poly(n) time; whereas for approximating P_f and F_f within an arbitrary multiplicative error, MCS cannot be run in poly(n) time.



 §.§ Main Ideas of the Proposed Algorithm

The proposed algorithm in Section <ref> provides an approximation of P_f and F_f in poly(N)-time with provable guarantees for any k-terminal reliability system. However, the number of cutsets (N) generally grows exponentially with the number of nodes (n). Thus, a natural question that arises is whether one can design an algorithm for computing (ϵ,δ)-approximations of P_f and F_f which runs in poly(n)-time. Karger in <cit.> proposed the first (and only) poly(n)-time algorithm for approximating P_f for all-terminal reliability systems (i.e., k=n). In this section, we present the first poly(n)-time algorithm for approximating F_f for all-terminal reliability systems. 

As was previously shown in <cit.>, in all-terminal reliability systems, the number of minimum cutsets, i.e., those cutsets with minimum weight, is poly(n). Moreover, there are a poly(n) number of near-minimum cutsets with weight not greater than a given constant factor of the weight of the minimum cutsets (see Lemma <ref>), and such cutsets can all be enumerated in poly(n) time (see Lemma <ref>). Note that, for a poly(n) number of cutsets, there still exist an exponential number of terms in (<ref>) and (<ref>), and consequently, using (<ref>) and (<ref>) one can only provide a series of upper- and lower-bounds on P_f and F_f via applying the bounding technique <cit.>. This, however, does not give poly(n)-time (ϵ,δ)-approximations of P_f and F_f, for arbitrary ϵ>0 and 0<δ<1. Using MCS, also, one may require exponentially many simulation runs to compute such approximations of P_f and F_f (see Section <ref>). 

To tackle this problem, we use the ideas and algorithms from the previous section, with new bounds to achieve the desired approximation parameters. In particular, we propose a fast and accurate algorithm to approximate F_f, using the near-minimum cutsets of weight no greater than α≥ 1 times the minimum cutset weight, for a proper choice of α upper bounded by 3+o(1) (i.e., using N^(α)≤ n^2α=O(n^6) cutsets (by Lemma <ref>)). Note that for approximating P_f, the proper choice of α, as was shown in <cit.>, is upper bounded by 2.  



 §.§ Proposed Algorithm

For an all-terminal reliability system, the inputs of the proposed algorithm are the failure rates {λ_i}, the repair rates {μ_i}, and the approximation parameters ϵ>0 and 0<δ<1. The algorithm proceeds in steps as follows:


  0. Initialization:


  0.1 Find a min-cutset 𝒞^*, and compute w^*=w(𝒞^*) and p^*=p(𝒞^*);

  0.2 Compute μ_min, λ_max, w_max, μ, μ^*=∑_i∈𝒞^*μ_i, s^* = min{max{w^*/w_max,1},m}, and ρ=μ_mins^*-λ_max (m-s^*);


  1. If p^*>n^-4: 


  1.1 Compute F̃^MC_f as the output of MCS for S = ⌈ (μ(2+ϵ)log 8)/ (p^*ρϵ^2)⌉ and T = ⌈ 12log(1/δ)⌉;

  1.2 Return F̃^MC_f.
 

  2. If p^*≤ n^-4:


  2.1 Take ξ = (ϵ/2)(ρ/μ)

  2.2 Take γ = (w^*/log n) -2, and α =1+(2/γ)+(log((2(γ+2)(μ-s^*μ_min))/(ξγ (μ-μ^*))))/(γlog n)

  2.3 Enumerate all N^(α) α-min cutsets;

  2.4 Compute a (ξ,δ/2)-approximation P̃^(α) of P using the KLM estimator; 

  2.5 Compute a (ξ,δ/2)-approximation P̃^(α)_f of P_f using the KLM estimator; 

  2.6 Return F̃_f = (P̃^(α)_f-P̃^(α))μ. 


The running time of this algorithm is O(n^4m(m+n)(1/ϵ^2)log(1/δ)) and O(N^(α) m^3(1/ϵ^2)log(1/δ)) for p^*>n^-4 and p^*≤ n^-4, respectively. (The running time of each step of the algorithm is given in Section <ref>.)

Theoretically, N^(α) = O(n^6) for our choice of α (see Lemma <ref>). This result follows from a worst-case analysis. However, for many practical systems, e.g., Internet2 network (see Section <ref>), N^(α) is much smaller, e.g., O(n^2). The threshold n^-4 for p^* in the algorithm is chosen to have a matching running time O(n^8) in both cases of p^* for systems with N^(α) = O(n^2). Nevertheless, for any choice of the threshold, the running time of the proposed algorithm is significantly less than that of MCS, for sufficiently small p^* (depending on the threshold's choice). 

The details of the enumeration of N^(α) α-min cutsets (Step 2.3) and the computations of P̃^(α)_f and P̃^(α) (Steps 2.4 and 2.5) are as follows. 



  §.§.§ Enumeration of α-min Cutsets
 We enumerate α-min cutsets by using a randomized algorithm, referred to as the  recursive generalized contraction (RGC) algorithm, due to Karger and Stein <cit.>. (The non-recursive and recursive original contraction algorithms output min-cutsets <cit.>.) For simplicity, we explain the non-recursive version of this algorithm.

The (non-recursive) generalized contraction algorithm proceeds in rounds. In each round, one component, say i, is randomly chosen with probability of choosing component i equal to w_i/w, and the two end-nodes of the component i are merged, while maintaining all components from either of these nodes to other nodes. The algorithm continues this process until more than ⌈ 2α⌉ nodes remain, and terminates otherwise. Once terminated, the algorithm selects a cutset in the resulting (multi-) system at random, and returns this cutset. 

Run the RGC algorithm for n^2αlog n^2α+c times for an arbitrary c>0 (in O(n^2αlog^2 n) time), and denote by C^(α) the set of all N^(α) (distinct) output cutsets.



  §.§.§ Computation of P̃^(α)
 Define Z_i and P_Z(i) as in Section <ref>, except using only the N^(α) cutsets in C^(α), instead of all the N cutsets in C; and construct a DNF formula 
    Φ_N^(α)=Z_1∨…∨ Z_N^(α).
 Run the KLM estimator with inputs (Φ_N^(α),P_Z;ξ,δ/2), and denote by P̃^(α) the output. 





  §.§.§ Computation of P̃_f^(α)



Similar to computing P̃^(α), in order to compute P̃^(α)_f, construct a DNF formula 
    Φ_N^(α),f=X_1∨…∨ X_N^(α),
 where X_i is defined as in Section <ref>, except only for i∈ [N^(α)], and not all i∈ [N]. Define P_X(i) accordingly as before. Run the KLM estimator with inputs (Φ_N^(α),f,P_X;ξ,δ/2), and denote by P̃^(α)_f the output. 






 §.§ Theoretical Analysis



The output of the algorithm in Section <ref>, F̃^MC_f or F̃_f, is an (ϵ,δ)-approximation of F_f. 







For the case of p^*>n^-4, the algorithm resorts to the MCS algorithm. We show that F̃^MC_f is an (ϵ,δ)-approximation of F_f. To this end, we use a version of the Chernoff's bound as follows. (The proof of this result can be found in the appendix.)


Let I_1,…,I_T be T independent and identically distributed random variables such that 0≤ I_t≤ 1 for all t. Let I≜∑_t=1^T I_t. Then, for any ϵ>0, 

    {I≥ (1+ϵ)𝔼(I)}≤exp(-(ϵ^2/2+ϵ)𝔼(I)),
 and 

    {I≤ (1-ϵ)𝔼(I)}≤exp(-(ϵ^2/2)𝔼(I)).



By the definition of φ_s,t in Step 2 of MCS, it is easy to see that 0≤ρ≤φ_s,t≤μ for all s,t (by (<ref>) and (<ref>)). Let φ̂_s,t≜φ_s,t/μ. Note that φ̂_s,t are independent and identically distributed random variables such that 0≤φ̂_s,t≤ 1 for all s,t. Recall that φ_t = (∑_s=1^Sφ_s,t)/S (Step 5 of MCS). Let φ̂_t ≜φ_t/μ. By the result of Lemma <ref>, for any ϵ>0, we get {|φ̂_t-𝔼(φ̂_t)|≥ϵ𝔼(φ̂_t)}≤ 2exp(-(ϵ^2/(2+ϵ))S𝔼(φ̂_t)) for all t. Note that 𝔼(φ_t)=F_f (by (<ref>) and (<ref>)), and accordingly, 𝔼(φ̂_t) = F_f/μ. Moreover, F_f≥ P_fρ (by (<ref>)) and P_f≥ p^*. (The probability that all components in a cutset are unavailable, P_f, is lower bounded by the probability that all components in a given min-cutset are unavailable, p^*, i.e., P_f≥ p^*.) Thus, for all t,

    {|φ_t-F_f|≥ϵ F_f}≤ 2exp(-(ϵ^2/2+ϵ) (ρ/μ) S p^*)	,
 or equivalently, 

    {|φ_t-F_f|≥ϵ F_f}≤1/4,
 for the choice of S in the algorithm (Step 1.1). Since F̃^MC_f, the output of MCS, is the median of {φ_t}_t=1^T, by applying the result of Lemma <ref>, it can be shown that 

    {|F̃_f^MC-F_f|≥ϵ F_f}   ≤δ
 for the choice of T in the algorithm (Step 1.1). (Similar technique is used in the appendix as part of the proof of Lemma <ref>.) Thus, F̃^MC_f is an (ϵ,δ)-approximation of F_f (by (<ref>)). (Similarly, MCS can compute such an approximation of P_f in poly(n) time.)
 
























Now, consider the case of p^*≤ n^-4. Note that, in this case, P_f can be arbitrarily small (in n and m), and S can be arbitrarily large. Thus, MCS cannot compute an (ϵ,δ)-approximation of F_f in poly(n) time. (A similar negative result holds for computing such an approximation of P_f.)

Fix an arbitrary α≥ 1. The following results hold for α-min cutsets. (The proofs are given in the appendix.) 


The number of α-min cutsets, N^(α), is bounded from above by n^2α.
 


Running the RGC algorithm n^2αlog n^2α+c times, for any c>0, one can enumerate all α-min cutsets in O(n^2αlog^2 n) time, with probability at least 1-n^-c.	


By the choice of the algorithm (Step 2.2), 

    γ=w^*/log n-2,
 i.e., p^*=n^-2-γ. (Since p^*≤ n^-4, it holds that γ≥ 2.) The probability that all components in a given cutset are unavailable is upper bounded by p^*, and the probability that they are all unexposed is upper bounded by 1-s^*(μ_min/μ). (Note that p(𝒞)≤ p^* and |𝒞|≥ s^* for all cutsets 𝒞.) Thus, the probability that all components in a given α-min cutset are unavailable and unexposed is upper bounded by p^*(1-s^*(μ_min/μ)), or equivalently, n^-2-γ(1-s^*(μ_min/μ)). 

By Lemma <ref>, there are at most n^2α α-min cutsets. By applying union bound, the probability that all components in some α-min cutset are unavailable and unexposed is upper bounded by n^-2-γ+2α(1-s^*(μ_min/μ)). This result is generalizable for all cutsets of weight greater than α w^* as follows. (A related, yet weaker, result was shown in <cit.>.) 



Let {𝒞_1,…,𝒞_N-N^(α)} be the set of all cutsets of weight greater than α w^*. Assume that w(𝒞_1)≤…≤ w(𝒞_N-N^(α)). Let M ≜min{n^2α, N-N^(α)}. First, consider the cutsets 𝒞_1,…,𝒞_M. For any 1≤ j≤ M, the probability that all components in 𝒞_j are unavailable and unexposed, exp(-w(𝒞_j)) (1-∑_i∈𝒞_jμ_i/μ), is upper bounded by (p^*)^α(1-s^*(μ_min/μ)). By applying union bound, the probability that all components in 𝒞_j for some 1≤ j≤ n^2α are unavailable and unexposed is upper bounded by 

    n^2α (p^*)^α(1- s^*μ_min/μ)=n^-αγ(1- s^*μ_min/μ).


Next, consider the remainder of the cutsets 𝒞_M+1,…,𝒞_N-N^(α) (if any). For any β>0, the number of cutsets of weight less than or equal to β w^* is upper bounded by n^2β (by Lemma <ref>). Thus, w(𝒞_n^2β)≥β w^*. This gives w(𝒞_j)≥ (w^*log j)/(2log n) for all j, and subsequently, the probability that all components in cutset 𝒞_j are unavailable and unexposed is upper bounded by (p^*)^(log j)/(2log n)(1-s^*(μ_min/μ))=j^-1-γ/2(1-s^*(μ_min/μ)). Again by a union-bound analysis, the probability that all components in cutset 𝒞_j for some j>n^2α are unavailable and unexposed is upper bounded by 

    ∑_j>n^2α j^-1-γ/2(1- s^*μ_min/μ)≤2/γ n^-αγ(1- s^*μ_min/μ).


Putting (<ref>) and (<ref>) together, the probability that all components in some cutset 𝒞_j, j∈ [N-N^(α)], are unavailable and unexposed is upper bounded by 

    n^-αγ(1-s^*μ_min/μ)+ 2/γ n^-αγ(1-s^*μ_min/μ) = n^-αγ(1+2/γ)(1-s^*μ_min/μ).
 This immediately yields the following result.  


The probability that all components in some cutset of weight greater than α w^* are unavailable and unexposed is bounded from above by 

    n^-αγ(1+2/γ)(1-s^*μ_min/μ).







The probability that all components in some cutset are unavailable and unexposed, P, is lower bounded by the probability that all components in the min-cutset 𝒞^*, found in Step 0.1, are unavailable and unexposed, p^*(1-μ^*/μ). Note that this argument holds for any cutset 𝒞, but the lower bound needs to be replaced with p(𝒞)(1-∑_i∈𝒞μ_i/μ). Thus, 

    P≥ p^*(1-μ^*/μ).
 (Note that the bound in (<ref>) can be improved as follows. Enumerate all N^(1)≤ n^2 min-cutsets, e.g., by using the RGC algorithm in O(n^2log^2 n) time, and select a min-cut 𝒞 with minimum ∑_i∈𝒞μ_i, e.g., via sorting in O(n^2log n) time.) 

By combining (<ref>) and (<ref>), it is easy to see that 

    n^-αγ(1+2/γ)(1-s^*μ_min/μ)≤ξ/2 P
 for any ξ>0 so long as 

    α≥1/γ(γ+2+log(2/ξ(γ+2/γ)(μ-s^*μ_min/μ - μ^*))/log n).
 By the choice of α in the algorithm (Step 2.2), it follows that (<ref>), and consequently, (<ref>) hold. 





Let P^(α) be the probability that all components in some α-min cutset are unavailable and unexposed. Then, 

    (1-ξ/2)P (a)≤ P^(α)(b)≤ P,
 where (a) follows from (<ref>), and (b) follows from the definitions of P^(α) and P. Thus, P^(α) is a (ξ/2,0)-approximation of P. Note that P^(α) corresponds to the α-min cutsets, and there are poly(n) such cutsets (Lemma <ref>), and they can be enumerated in poly(n) time (Lemma <ref>). However, one cannot compute P^(α) in poly(n) time via the inclusion-exclusion formula due to the exponential number of terms. We, instead, approximate P^(α) in poly(n) time using the KLM estimator. 

The formula Φ_N^(α) is true so long as some clause Z_i is true. Thus, the truth probability of Φ_N^(α) is equal to P^(α). Since Φ_N^(α) is a DNF formula, the output of the KLM estimator with inputs (Φ_N^(α),P_Z;ξ/2,δ/2), denoted by P̃^(α), is a (ξ/2,δ/2)-approximation of P^(α), i.e., 

    {|P̃^(α)-P^(α)|≥ξ/2 P^(α)}≤δ/2.
 By using (<ref>) and (<ref>), we get  

    {|P̃^(α)-P|≥ξ P }≤δ/2.
 This immediately yields the following result.


P̃^(α) is a (ξ,δ/2)-approximation of P.


Let P_f^(α) be the probability that all components in some α-min cutset are unavailable. Similar to (<ref>), as was shown in <cit.>, it follows that

    (1-ξ/2)P_f ≤ P_f^(α)≤ P_f.
 

The probability that the formula Φ_N^(α),f is true is equal to P^(α)_f. Thus, the output of the KLM estimator with inputs (Φ_N^(α),f,P_X;ξ/2,δ/2), denoted by P̃_f^(α), is a (ξ/2,δ/2)-approximation of P^(α)_f. Thus, 

    {|P̃_f^(α)-P_f^(α)|≥ξ/2 P_f^(α)}≤δ/2.
 By using (<ref>) and (<ref>), we get

    {|P̃_f^(α)-P_f|≥ξ P_f }≤δ/2.
 Then, the following result is immediate. 


P̃_f^(α) is a (ξ,δ/2)-approximation of P_f.


Putting together (<ref>) and (<ref>), it follows that

    {|(P̃_f^(α)-P̃^(α))-(P_f-P)|≥ξ (P_f+P) }≤δ.
 The rest of the proof is the same as that in the proof of Theorem <ref> (and hence omitted), except that P̃_f and P̃ are replaced with P̃^(α)_f and P̃^(α).

Similarly as in (<ref>), it can be shown that

    {|F̃_f^MC-F_f|≥ϵ}   ≤δ
 so long as 

    S≥μ(2F_f+ϵ/ϵ^2)log 8
 and 

    T≥ 12log(1/δ).

Taking S = ⌈ (μ(2μ+ϵ)log 8)/ϵ^2⌉ (noting F_f≤ P_fμ≤μ (by (<ref>))) and T = ⌈12log(1/δ)⌉, MCS can approximate F_f, with an additive error of at most ϵ and an error probability of at most δ, in poly(n) time (by (<ref>)-(<ref>)). (Since μ=O(m)=O(n^2), then ST=O((n^4/ϵ^2+n^2/ϵ)log(1/δ)).) Similar result holds for approximating P_f using MCS. Note, however, that MCS cannot compute approximations of P_f and F_f within a multiplicative error in poly(n) time. This suggests that approximating P_f and F_f within a multiplicative error is more challenging than that within an additive error. 





 §.§ Computational Complexity

The initialization (Step 0) can be run in O(n^3) time. Step 0.1 can be run in O(n^3) time, e.g., by finding a min-cutset in O(nm+n^2log n) time by the Stoer-Wagner algorithm <cit.>, and computing p^* and w^* in O(m) time; In Step 0.2, μ_min, λ_max, and w_max can be computed in O(mlog m) time (via sorting), μ and μ^* in O(m) time, and s^* and ρ in O(1) time. 



For simplicity, we ignore ϵ and δ in what follows, since they do not depend on n and m. In Step 1, F̃_f^MC can be computed in O(n^4m(m+n))=O(n^8) time by using MCS. Step 2 can also be run in O(n^12) time as follows. (These running times follow from a worst-case analysis, and some of the bounds are far from tight.) Since γ≥ 2 and ξ=Ω(n^-2) for our choice of ξ (Step 2.1), it is easy to see that our choice of α (Step 2.2) is upper bounded by 3+O(1/log n), or equivalently, 3+o(1). Moreover, N^(α)≤ n^2α (by Lemma <ref>). Thus, in Step 2.3, it suffices to enumerate N^(α)=O(n^6) α-min cutsets for the choice of α in the algorithm. (For computing P̃_f^(α), as was shown in <cit.>, it suffices to enumerate α-min cutsets for some α≤ 2, and there are O(n^4) such cutsets.) Step 2.3 can be run, e.g., using the RGC algorithm <cit.>, in O(n^2αlog^2 n)=O(n^6log^2 n) time (by Lemma <ref>). In Steps 2.4 and 2.5, the estimates P̃^(α) and P̃^(α)_f can be computed by running the (self-adjusting) KLM estimator in O(N^(α) m^3) = O(n^12) time. From these arguments, we conclude that the proposed algorithm runs in poly(n) time.



§ SIMULATION RESULTS

In this section, we compare the accuracy and the running time of the Monte Carlo simulation and the proposed algorithm for approximating F_f for a 3× 3 grid network (Fig. <ref>) and for the Internet2 network <cit.>, particularly, the layer 3 (Fig. <ref>) and the layer 2 (Fig. <ref>). We refer to these systems as S1, S2, and S3, respectively. S1, S2, and S3 contain n=9, 20, and 35 nodes and m=12, 31, and 46 components, respectively. 



For simulations, we have used the Matlab Parallel Computing toolbox for 32 workers on CentOS linux machines with 2X - Intel Xeon E5-2697A V4 @ 2.6GHz processors and 512GB memory. 





The unavailability probability of all components in the system S1 is assumed to be p for different values of p∈{10^-2,10^-2.2,…,10^-3.8}. The unavailability probability of each link (component) in the systems S2 and S3 is assumed to be proportional to the number of hops between the two end-points of the link in the layer 1 (see, for details, <cit.>). For example, for a given p, the unavailability probability of the link between Seattle and Salt Lake City in the system S2 is assumed to be 2p since they are two hops away in the layer 1, and the unavailability probability of the link between Denver and El Paso in the system S3 is assumed to be 3p since they are three hops away in the layer 1. For the systems S2 and S3, we consider different values of p∈{10^-3,10^-3.2,…,10^-4.8} and p∈{10^-4,10^-4.2,…,10^-5.8}, respectively. Note that we consider cases with very small values of p≪ 1 since the failure of the system is a rare event in such cases. 

For each component with unavailability probability θ p (for some integer θ), the repair rate and the failure rate of the component are assumed to be 1 and θ p/(1-θ p), respectively. 





Table <ref> lists the maximum failure probability of a cutset (p^*), the parameter α, and the number of α-min cutsets being enumerated (N^(α)), for each of the systems S1, S2, and S3 and each p. To enumerate the α-min cutsets, for each system and each p, we ran the RGC algorithm 2× 10^4 times, and the total running time was about 30 seconds for each case. 






For fixed δ = 10^-2 (i.e., an error probability at most 1%), Table <ref> represents the approximations F̃_f and F̃^MC_f of F_f for the system S1 using the proposed algorithm and the MCS algorithm, the total running time for each algorithm, and the theoretical error factor and the actual (observed) error factor of each algorithm. In particular, the theoretical error factor (ϵ) is computed based on the analysis in Section <ref>, and the actual error factor, i.e., max{|φ-F^-_f|/F^-_f,|φ-F^+_f|/F^-_f}, where φ is F̃_f or F̃^MC_f for the proposed or the MCS algorithm, respectively, is computed based on the first-order upper- and lower-bounds F_f^+ and F_f^- on F_f resulting from the bounding technique. The values of F_f^+ and F_f^- are also given for reference in Table <ref>. Similarly, Tables <ref> and <ref> correspond to the results for the systems S2 and S3, respectively.

To compute F_f^+ and F_f^-, we enumerated all cutsets in each system. The accuracy of this technique, however, cannot be fairly compared with that of the other two algorithms since the running time for enumerating all cutsets in each system (about 1 week for the system S2 and 2 weeks for the system S3) was much larger than that of the other two algorithms. 

A simple comparison of the MCS algorithm and the proposed algorithm shows that not only is the latter technique faster than the former, but also it provides more accuracy, for each system and for any given p. For example, for the system S1 and p=10^-3, running the proposed algorithm for 181 seconds one can approximate F_f with an error factor of 2.81× 10^-5; whereas running the MCS algorithm for 289 seconds, one can only approximate F_f with an error factor of 1.05× 10^-2. As an another example, for the system S1 and p=10^-3.8, the proposed algorithm provides an approximation with an error factor of 4.65× 10^-7 in 214 seconds; whereas running the MCS algorithm for 333 seconds, no failure state is detected (and hence the output is “zero,” and the actual error factor is meaningless). Similar comparison results hold for larger systems S2 and S3 (see Tables <ref> and <ref>). Note that the advantages of the proposed algorithm over the MCS algorithm are even more evident for larger systems.

These results confirm that the proposed algorithm offers a significantly better tradeoff between running time and approximation accuracy. Note also that the improvements become more profound for smaller values of p. This is evident from the two examples above. Furthermore, for smaller values of p the running time of the proposed algorithm becomes smaller, whereas the running time of the MCS algorithm remains almost the same. This comes from the fact that as p decreases the parameter α decreases (even for smaller error factor), and the number of α-min cutsets decreases. 

Furthermore, the theoretical error factors for both algorithms, obtained from the worst-case analysis in Section <ref>, are much larger than the actual error factors, computed based on the simulation results. This suggests that for many practical systems each of these algorithms may provide a better complexity-accuracy tradeoff than what was shown in the analysis. However, the ratio of the actual error factor to the theoretical error factor for the proposed algorithm, when compared to that for the MCS algorithm, is much smaller. Thus, one can expect that the proposed algorithm outperforms the MCS algorithm in practice even more than that in theory.   




§ CONCLUSION AND OPEN PROBLEMS

In this work, we considered the problem of estimating the failure frequency of large-scale composite k-terminal reliability systems. It was previously shown that the failure probability can be efficiently approximated with provable guarantees. However, no such result was previously known for approximating the failure frequency. 

We proposed the first polynomial-time (in the number of cutsets in the system) algorithm for approximating the failure frequency within an arbitrary multiplicative error and with an arbitrary error probability. The main ideas of the proposed algorithm are in summary as follows: (i) the failure frequency can be linked to a linear combination of the probabilities of two sets of events, and (ii) each of these probabilities can be written as a Boolean formula of disjunctive normal form, and the truth probability of such a formula can be estimated using an unbiased estimator within an arbitrary error factor and with an arbitrary error probability. 

The number of cutsets of a system can generally grow exponentially with the number of nodes in the system. Motivated by this, for the special case of all-terminal reliability systems in which all nodes are terminals, we proposed the first polynomial-time (in the number of nodes in the system) algorithm using only the near-minimum cutsets, instead of all cutsets, for approximating the failure frequency with an arbitrary multiplicative error factor and with an arbitrary error probability. The main ideas here are: (i) the number of near-minimum cutsets, as opposed to the number of all cutsets, is only polynomial in the number of nodes, and all such cutsets can be enumerated in polynomial time, and (ii) for a proper choice of near-minimum cutsets, neglecting all those cutsets which are not near-minimum yields an arbitrarily small error while approximating the failure frequency. 

We compared the accuracy and the running time of the Monte Carlo simulation (MCS) and the proposed algorithm for all-terminal reliability systems (via theoretical analysis and simulation study). The comparison results confirm that the proposed algorithm achieves a significantly better tradeoff between accuracy and running time than MCS.

Unlike the special case of all-terminal reliability systems, the number of near-minimum cutsets in k-terminal reliability systems (for arbitrary 2≤ k<n) generally grows exponentially with the number of nodes. One, and perhaps the most important, problem which remains open in this area of research is to design a polynomial-time (in the number of nodes) algorithm for approximating the failure probability and the failure frequency in k-terminal reliability systems. Some other directions for future research include deriving tighter bounds for the analysis and improving the running time of the proposed algorithms. 





 §.§ Proof of Lemma <ref>

Let P(z) be the probability that a randomly chosen assignment is equal to z. Given that Z_j is chosen in Step 1, the probability of choosing z in Step 2 is equal to P(z)/P_Z(j). By the linearity of expectation, it follows that 

    𝔼(π_s,t| Z_j)=∑_z: z Z_jQ_Z/N(z)P(z)/P_Z(j),
 where the notation “z Z_j” indicates that the summation is taken over all z satisfying Z_j. (Similarly, in the following, we use the notation “zΦ_M” to indicate all z satisfying Φ_M.) This gives 

    𝔼(π_s,t)    (a)=   ∑_jP_Z(j)/Q_Z∑_z: z Z_jQ_Z/N(z)P(z)/P_Z(j)
       (b)=   ∑_j∑_z: z Z_jP(z)/N(z)
       (c)=   ∑_z: zΦ_M∑_j: z Z_jP(z)/N(z)
       (d)=   ∑_z: zΦ_MP(z)/N(z)∑_j: z Z_j 1
       (e)=   ∑_z: zΦ_M P(z)
       (f)=   Π_M,
 where (a) follows from the law of total expectation, i.e., 
    𝔼(π_s,t)=∑_j{selecting  Z_j}·𝔼(π_s,t| Z_j),
 noting that the probability of selecting Z_j is equal to P_Z(j)/Q_Z, and 𝔼(π_s,t|Z_j) is given by (<ref>); (b) follows since Q_Z and P_Z(j) are independent of z; (c) follows since any z satisfying Z_j (for any j) also satisfies Φ_M; (d) follows since P(z)/N(z) is independent of j; (e) follows since ∑_j: z Z_j 1 = N(z) (by definition), and (f) follows from the definition of Π_M. Thus, 
    𝔼(π_s,t)=Π_M.
 Similarly, 

    𝔼((π_s,t)^2)    =   ∑_jP_Z(j)/Q_Z∑_z: z Z_j(Q_Z)^2/(N(z))^2P(z)/P_Z(j)
       =   ∑_j∑_z: z Z_jQ_Z/(N(z))^2 P(z)
       =   ∑_z: zΦ_MQ_Z/N(z) P(z)
       ≤    Q_Z Π_M,
 noting that N(z)≥ 1 for all z satisfying Φ_M. Thus,  

    var(π_s,t)	    =   𝔼((π_s,t)^2)-(𝔼(π_s,t))^2
       ≤    Q_Z Π_M - (Π_M)^2
       ≤    (M-1)(Π_M)^2,
 noting that 

    Q_Z    =∑_j P_Z(j) 
       = ∑_j{a randomly chosen z satisfies  Z_j}
       ≤  M ·{a randomly chosen z satisfies Φ_M}
       = MΠ_M.
 Since π_t = (∑_sπ_s,t)/S for all t∈ [T], it follows that 𝔼(π_t) = 𝔼(π_s,t) and var(π_t)=var(π_s,t)/S. Thus,

    𝔼(π_t) =Π_M
 and 
    var(π_t)≤(M-1) (Π_M)^2/S
 for any t∈ [T]. Applying the Chebychev's inequality, we get

    {|π_t -Π_M| ≥ξΠ_M}
       ≤{|π_t-Π_M|≥ξ√(S/M-1var(π_t) )}
       ≤M-1/ξ^2 S
       ≤1/4,
 for the choice of S in the algorithm (Step 5). Define an indicator variable I_t for all t∈ [T] as follows: I_t=1 if |π_t-Π_M|≥ξΠ_M, and I_t=0 otherwise. Let I ≜∑_t=1^T I_t. Note that 𝔼(I)≤ T/4 since I_t=1 with probability at most 1/4, and I_t=0 otherwise. By applying Lemma <ref>, we get 
    {I≥ (1+ϵ)𝔼(I)}≤exp(-(ϵ^2/2+ϵ)𝔼(I))
 for any ϵ>0. Taking 
    ϵ = T-2𝔼(I)/2𝔼(I),
 we have (1+ϵ)𝔼(I) = T/2. Since 𝔼(I)≤ T/4, then T-2𝔼(I)≥ T/2 and T+2𝔼(I)≤ 3T/2. Thus, 

    {I≥ (1+ϵ)𝔼(I)}
       = {I≥T/2}
       ≤exp(-(ϵ^2/2+ϵ)𝔼(I))
       = exp(-1/2((T-2𝔼(I))^2/T+2𝔼(I)))
       ≤exp(-T/12)
       ≤δ,
 for the choice of T in the algorithm (Step 7). Note that Π̃_M=median({π_t}_t=1^T) (by definition). If the median Π̃_M deviates from Π_M by more than ξΠ_M (i.e., |Π̃_M-Π_M|≥ξΠ_M), then π_t deviates from Π_M by more than ξΠ_M (i.e., |π_t-Π_M|≥ξΠ_M) for at least half of t∈ [T], or equivalently, I_t=1 for at least half of t∈ [T] (i.e., I≥ T/2). Thus, 

    {|Π̃_M-Π_M|≥ξΠ_M}
       ≤{I ≥T/2}.
 By combining (<ref>) and (<ref>), it is easy to see that Π̃_M is a (ξ,δ)-approximation of Π_M.









 §.§ Proof of Lemma <ref>

The following inequality is useful for the proof of the lemma. For any x> 0 and 0≤θ≤ 1,

    x^θ≤ 1+(x-1)θ.
 This inequality follows immediatley from a generalization of the Brnouolli's inequality as follows. 

For any x>-1 and 0≤θ≤ 1, 
    (1+x)^θ≤ 1+xθ.

 

Taking x = e^z for any z, we get e^z I_t≤ 1+(e^z-1)I_t for all t (by (<ref>)). Taking expectation from both sides of this inequality, it follows that 𝔼(e^z I_t)≤𝔼(1+(e^z-1)I_t) = 1+(e^z-1)𝔼(I_t). Note that 𝔼(I_t) = 𝔼(I)/T since I = ∑_t=1^T I_t, and I_1,…,I_T are identically distributed. Thus, one can see that

    𝔼(e^z I_t)≤ 1+(e^z-1)𝔼(I)/T.
 Similarly, it can be seen that

    𝔼(e^-z I_t)≤ 1+(e^-z-1)𝔼(I)/T.
 Obviously, 

    {I≥ (1+ϵ)𝔼(I)} = {e^z I≥ e^z(1+ϵ)𝔼(I)},
 for any z>0. By the Markov's inequality, {e^z I≥ e^z(1+ϵ)𝔼(I)}≤𝔼(e^z I)/e^z(1+ϵ)𝔼(I). Since I_1,…,I_T are independent, 𝔼(e^z I) = 𝔼(e^z I_1)⋯𝔼(e^z I_T). By (<ref>), 𝔼(e^z I)≤ (1+(e^z-1)𝔼(I)/T)^T. Since (1+x/T)^T≤ e^x for any x, then 𝔼(e^z I)≤ e^(e^z-1)𝔼(I). Thus, 

    {e^z I≥ e^z(1+ϵ)𝔼(I)}   ≤e^(e^z-1)𝔼(I)/e^z(1+ϵ)𝔼(I)
        = e^(e^z-1-z(1+ϵ))𝔼(I).
 Let f(z)≜ e^z-1-z(1+ϵ). Taking z = log(1+ϵ), we get f(z)=(1+ϵ)(1-log(1+ϵ))-1. It is easy to see that log(1+ϵ)≥ 2ϵ/(2+ϵ) for any ϵ≥ 0. (Letting g(ϵ)≜log(1+ϵ)-2ϵ/(2+ϵ), and noting g(0)=0 and g(ϵ) is increasing, it follows that g(ϵ)≥ 0 for any ϵ≥ 0.) Thus, 

    f(z)≤ -ϵ^2/(2+ϵ)

for any ϵ≥ 0. By (<ref>)-(<ref>), we get 

    {I≥ (1+ϵ)𝔼(I)}≤ e^-(ϵ^2/2+ϵ)𝔼(I).
 Similarly as in (<ref>), it is obvious that

    {I≤ (1-ϵ)𝔼(I)} = {e^-z I≥ e^-z(1-ϵ)𝔼(I)},
 for any z>0. Similar to (<ref>), except by using (<ref>) instead of (<ref>), it follows that 
 
    {e^-z I≥ e^-z(1-ϵ)𝔼(I)}≤ e^(e^-z-1+z(1-ϵ))𝔼(I).
 Let f(z)≜ e^-z-1+z(1-ϵ). First, suppose that ϵ≥ 1. Taking z=ϵ, we get f(z)=e^-ϵ-1+ϵ(1-ϵ). Since e^-ϵ≤ 1/(1+ϵ) and ϵ/(1+ϵ)≥ 1/2, then f(z)≤ -ϵ^2(ϵ/(1+ϵ))≤ -ϵ^2/2 for any ϵ≥ 1. Next, suppose that 0≤ϵ<1. Taking z = -log(1-ϵ), we get f(z)=-(1-ϵ)log(1-ϵ)-ϵ. By the Taylor's expansion, log(1-ϵ)=-ϵ-ϵ^2/2-ϵ^3/3-ϵ^4/4-…, and consequently, (1-ϵ)log(1-ϵ)=-ϵ+ϵ^2/2+ϵ^3/6+ϵ^4/12+…≥ -ϵ+ϵ^2/2. Then, f(z)≤ϵ-ϵ^2/2-ϵ= -ϵ^2/2 for any 0≤ϵ<1. By these arguments, 
    f(z)≤ -ϵ^2/2

for any ϵ≥ 0. By (<ref>)-(<ref>), we get 

    {I≤ (1-ϵ)𝔼(I)}≤ e^-(ϵ^2/2)𝔼(I).
 

   
  



 §.§ Proof of Lemma <ref>

Fix an arbitrary α≥ 1. Consider an application of the (non-recursive) generalized contraction algorithm. Recall that in each round of this algorithm, one component, say i, is chosen at random with probability of choosing component i equal to w_i/w, and the two end-nodes of the component i are merged. The algorithm continues this process until more than ⌈ 2α⌉ nodes remain, and terminates otherwise. Once terminated, the algorithm returns a randomly selected cutset in the resulting (multi-) system. 

Assume, without loss of generality, that the weights of all components are the same and equal to w_0. Otherwise, we can replace each component i of weight w_i by w_i/w_0 parallel components, each of weight w_0, for sufficiently small w_0 such that w_i/w_0 is an integer for all i (for more details, see <cit.>). Note that the minimum weight of a cutset in this new system is equal to that in the original system, whereas the minimum size of a cutset in this new system can be different from that in the original system. We denote by w^* and s^* (=w^*/w_0), with a slight abuse of notation, the minimum weight and the minimum size of a cutset in the new system, respectively. 



Let 𝒞 be an arbitrary cutset of weight β w^* for some 1≤β≤α. Note that 𝒞 has β s^* components. We say that 𝒞 is hit in round i if one of its components is chosen and collapsed in round i. Consider the round n-r+1 where r nodes are remaining (for arbitrary ⌈ 2α⌉< r≤ n). Since the number of components connected to each node is lower bounded by s^* (otherwise there exists a cutset of size less than s^*), then the total number of components is lower bounded by r s^*/2 (otherwise there exists a node with less than s^* components connected to it). Thus the probability that 𝒞 is hit in round n-r+1 is upper bounded by β s^*/(r s^*/2)=2β/r≤ 2α/r. Similarly, the probability that 𝒞 survives (i.e., 𝒞 is not hit in) rounds 0,1,…,i, and 𝒞 is hit in round i+1 is upper bounded by β s^*/((n-i) s^*/2)=2β/(n-i) for any 0≤ i≤ n-⌈ 2α⌉-1. Thus, the probability that 𝒞 survives all rounds until more than ⌈ 2α⌉ nodes remain is lower bounded by 

    ∏_i=0^n-⌈2α⌉ -1(1-2α/n-i) < (2/n)^2α.
 (See, for more details, the proof of <cit.>.) 

Once the algorithm terminates, there exist ⌈ 2α⌉ nodes in the eystem, and consequently, the number of cutsets is upper bounded by the number of bipartitions of the nodes, i.e., 2^⌈ 2α⌉-1-1<2^2α. Thus, the probability that 𝒞 is chosen via the random selection is lower bounded by 1/2^2α. 

The probability that 𝒞 is output by the algorithm is the product of two probabilities: (i) the probability that 𝒞 survives until ⌈ 2α⌉ nodes remain (this probability is lower bounded by 2^2α/n^2α), and (ii) the probability that 𝒞 is chosen by the random selection (this probability is lower bounded by 1/2^2α). Thus the probability that 𝒞 is output by the algorithm is lower bounded by (2^2α/n^2α) × (1/2^2α)= 1/n^2α. 

Each run of the algorithm outputs an α-min cutset with probability at least 1/n^2α. Thus, it follows that the number of α-min cutsets is at most n^2α. (This holds because the algorithm returns each distinct α-min cutset with probability lower bounded by 1/n^2α, and such events are disjoint.) 

Thus O(n^2) runs of the contraction algorithm are sufficient to find a min-cutset with high probability. By a clever recursive implementation of the contraction algorithm, as shown in <cit.> and <cit.>, a min-cutset can be found in O(n^2log^2 n) time (instead of O(n^4) time for the obvious implementation) with high probability.



 §.§ Proof of Lemma <ref>

The proof follows from the coupon-collector argument <cit.>: if there are M bins, and potentially an infinite number of balls to be thrown independently and uniformly one at a time, then throwing Mlog (M/δ) balls suffices with probability at least 1-δ to have each bin contain at least one ball. (To be specific, throwing Mlog M - Mloglog (1/(1-δ)) balls suffices with probability 1-δ as M grows large <cit.>.) By the result of Lemma <ref>, the number of α-min cutsets is at most n^2α. Think of α-min cutsets as bins and the runs of the (non-recursive) generalized contraction algorithm as balls. By the coupon-collector argument, one can enumerate all α-min cutsets with probability at least 1-n^-c, for any c>0, in O(n^2α+2log(n^2α+c)) time  (by running the generalized contraction algorithm n^2αlog (n^2α+c) times). By running RGC algorithm n^2αlog (n^2α+c) times, as was shown in <cit.> and <cit.>, all α-min cutsets can be found in O(n^2α(log (n^2α+c))(log n)) = O(n^2αlog^2 n) time with probability at least 1-n^-c (see, e.g., <cit.>). 




IEEEtran


