Efficient presolving methods for influence maximization problem in social networks
Sheng-Jie Chen, Wei-Kun Chen, Yu-Hong Dai, Jian-Hua Yuan, Hou-Shan Zhang
EEfficient presolving methods for influence maximizationproblem in social networks
Sheng-Jie Chen , Wei-Kun Chen ,Yu-Hong Dai , Jian-Hua Yuan , Hou-Shan Zhang January 2, 2021
Abstract
The influence maximization problem (IMP) is a hot topic, which asks foridentifying a limited number of key individuals to spread influence in a socialnetwork such that the expected number of influenced individuals is maximized. Thestochastic maximal covering location problem (SMCLP) formulation is a mixedinteger programming formulation that approximates the IMP by the Monte-Carlosampling. However, the SMCLP formulation cannot be solved efficiently usingexisting exact algorithms due to its large problem size. In this paper, we concentrateon deriving presolving methods to reduce the problem size and hence enhance thecapability of employing exact algorithms in solving the IMP. In particular, wepropose two effective presolving methods, called the strongly connected nodesaggregation (SCNA) and the isomorphic nodes aggregation (INA). The SCNAenables us to build a new SMCLP formulation that is much more compact thanthe existing one. For the INA, an analysis is given on the one-way bipartite socialnetwork. Specifically, we show that under certain conditions, the problem size of thereduced SMCLP formulation depends only on the size of the given network but noton the number of samplings and this reduced formulation is strongly polynomial-time solvable. Finally, we integrate the proposed presolving methods SCNA andINA into the Benders decomposition algorithm, which is recognized as one of thestate-of-the-art exact algorithms for solving the IMP. Numerical results demonstratethat with the SCNA and INA, the Benders decomposition algorithm is much moreeffective in solving the IMP in terms of solution time.
Keywords
Benders decomposition · Influence maximization · Integer program-ming · Presolving methods · Social network · Stochastic programming
Mathematics Subject Classification · Nowadays, with the popularity of online social networks such as Facebook, Instagram,and Twitter, propagation of influence in social networks has received more and moreattention. Promotion of products, ideas, and specific behavior patterns can all be viewedas propagation of influence. In practice, influence spreads among individuals through the Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China;School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China, { shengjie_ chen,dyh}@ lsec. cc. ac. cn School of Mathematics and Statistics/Beijing Key Laboratory on MCAACI, Beijing Institute ofTechnology, Beijing 100081, China, chenweikun@ bit. edu. cn School of Science, Beijing University of Posts and Telecommunications, Beijing 100876, China, { houshan_ zhang,jianhuayuan}@ bupt. edu. cn a r X i v : . [ m a t h . O C ] J a n o-called “word-of-mouth” exchanges. Individuals with more social connections can beseen as more influential, which means they are more likely to exert influence on others. Inthis setting, one related optimization problem, called the influence maximization problem (IMP), is to select a limited number of key individuals as seed set S to trigger a spreadprocess in the social network such that the expected number of influenced individuals ismaximized after the spread. Mathematically, the IMP can be written asmax S⊆V , | S |≤ K µ ( S ) , (1)where V is the set of individuals in the social network, S ⊆ V ( |S| ≤ K ) is the seed set ofkey individuals that need to be identified, and µ ( S ) is the influence function measuringthe expected number of individuals in the social network that can be influenced by theindividuals in seed set S . Kempe et al. [17] first proposed the discrete optimization problem formulation (1)for the IMP. Depending on different diffusion ways, they introduced two fundamentalpropagation models for the IMP: the independent cascade model (ICM) and the linearthreshold model (LTM). They showed that under both the ICM and LTM, the IMPis NP-hard, indicating that obtaining an optimal solution is challenging for the large-scale cases. However, by proving that the influence function µ ( S ) is monotone andsubmodular, they were able to design a greedy algorithm, which starts with S = ∅ and iteratively adds the individual with maximal marginal gain, with an approximationratio of (1 − /e ) where e is the base of the natural logarithm. Unfortunately, givena fixed seed set S , it is hard to compute µ ( S ) exactly (indeed, as shown in [4, 5], thisproblem is P-hard). Therefore, the Monte-Carlo sampling approach is used in [17]to approximate the influence spread and estimate influence function µ ( S ). Following[17], many researchers focused on improving the greedy algorithm and developing otherheuristic algorithms for solving the IMP. Specifically, Leskovec et al. [20] utilized thesubmodularity of µ ( S ) and presented an improved greedy algorithm called cost-effectivelazy forward . According to their numerical results, their method is almost 700 timesfaster than the basic greedy algorithm of Kempe et al. [17]. Chen et al. [3] proposedanother algorithmic enhancement for the greedy algorithm which reduces the graphsearching time on computing the marginal increment of a given individual. Furthermore,they developed a much more efficient algorithm, called degree discount heuristic, for theIMP under the ICM that nearly matches the performance of the greedy algorithm. The two-phase influence maximization [32] and the influence maximization via martingales [33] heuristic algorithms also deserve special attention. They can not only guaranteean approximation ratio of (1 − /e ), but also enable to solve large-scale IMPs in nearlylinear time. We refer to [18, 4, 6, 10] for more greedy or heuristic algorithms for solvingthe IMP and [24] for a detailed comparison among different heuristic algorithms.Most of the greedy and other heuristic algorithms find a suboptimal solution for theIMP with some worst case guarantees. However, in some applications, it is very crucialto identify an optimal solution instead of just a suboptimal one; see [11]. As a result,using exact algorithms to solve the IMP has attracted more and more attentions recently.In most of the exact algorithms, a mixed integer programming (MIP) formulation isestablished. In particular, Wu et al. [34] transferred the IMP into the so-called two-stage stochastic submodular MIP model , and proposed a delayed constraint generation algorithm to solve the problem to optimality. The computational results indicate thattheir algorithm is more efficient than the basic greedy algorithm in [17], especially when K is large. G¨uney [11] and Li et al. [23] formulated the IMP as a stochastic maximalcovering location problem (SMCLP) which is built on the live-arc graphs where each of2hem corresponds to a scenario. The SMCLP can be solved to global optimality usingstandard MIP solvers like CPLEX [7] and GUROBI [14]. Along this line, G¨uney et al.[12] developed a reformulation of the SMCLP and proposed a Benders decomposition (BD) algorithm. Their experiment results show that the BD algorithm outperforms theone in [34] by several orders of magnitude in terms of solution time. For more exactalgorithms of the IMP and its variants, we refer to [9, 28, 13, 16].However, due to the NP-hardness of the IMP , the above exact algorithms are stillinefficient, especially when the size of the social network is large.
Presolving [1] isan appealing strategy to address this issue. It removes redundant information andstrengthens the model formulation with the aim of improving the performance of thesubsequent solution procedure (e.g., the branch-and-cut or the BD approach). Indeed,presolving has been recognized as a standard routine of the state-of-the-art MIP solvers(e.g., CPLEX [7] and GUROBI [14]). For the problems with specific structures, developingcustomized presolving methods is often much more effective; see [2, 15, 25] for usingvarious customized presolving methods to solve some specific problems. In terms of theIMP, few articles are devoted to designing customized presolving methods. To the best ofour knowledge, only two simple presolving methods have been developed in the literature[12, 16] and they have been proved to be beneficial to solving the IMP in certain cases.Consequently, developing more customized presolving methods to further enhance thecapability of using exact algorithms to solve the IMP is still highly needed.
In this paper, we attempt to develop more presolving methods for the IMP based on theSMCLP formulation. The main contributions of this paper are summarized as follows. • By exploiting the problem structure of the IMP, we propose two new presolvingmethods including (i) the strongly connected nodes aggregation (SCNA) which aggre-gates nodes in each strongly connected component (SCC), in the same live-arc graph,into a single (virtual) node; and (ii) the isomorphic nodes aggregation (INA), whichextends the above idea into using isomorphic nodes among different live-arc graphsfor aggregation (two nodes in different live-arc graphs are isomorphic if the nodesthat can influence them are the same in the associated live-arc graphs). • We provide a theoretical analysis of the SCNA and INA to illustrate their strengthsto reduce the problem size of the IMP. In particular, for the SCNA, we show thatafter the exhaustive applications of this method, the SMCLP formulation can be builton the compact live-arc graphs (obtained by aggregating all SCCs in the originallive-arc graphs) and hence leads to a much smaller problem formulation. As for theINA, we consider a special case that the IMP is built based on a one-way bipartitenetwork and under the LTM. In this case, we show that after presolving, the numbersof variables and constraints are linear with the size of the network (independent ofthe number of samplings). Furthermore, the linear programming (LP) relaxation isshown to be tight and the problem is proved to be strongly polynomial-time solvable. • We integrate the SCNA and INA into the BD algorithm [12], which is recognizedas one of the state-of-the-art exact algorithms to solve the IMP. We show that theproposed SCNA and INA enable us to develop a much faster separation algorithmfor the Benders cuts, as compared with the one in [12].Our extensive numerical results on real-world social networks demonstrate that (i) theproposed SCNA and INA are indeed quite effective in reducing the problem size ofthe SMCLP formulation; (ii) when integrating them into the BD algorithm, they cansignificantly speed up the solution procedure of the IMP.The remainder of the paper is organized as follows. Section 2 briefly reviews the twomain propagation models (ICM and LTM) and the SMCLP formulation for the IMP.3ection 3 presents the two presolving methods and some theoretical analysis. Section 4describes the implementation of the proposed presolving methods and the integration ofthem with the BD algorithm for solving the IMP. Section 5 provides the computationalresults. Finally, section 6 gives some concluding remarks.
In this section, we briefly review the propagation models and the SMCLP formulation forthe IMP. We use a directed graph G = ( V , A ) to refer to a social network, in which a node i ∈ V represents the individual involved in the influence spread, and an arc ( i, j ) ∈ A represents individual i has the potential ability to influence (or activate) individual j .The spread of influence in a given social network G needs to obey certain propagationrules. In [17], Kempe et al. provided the following two fundamental influence spreadmodels called ICM and LTM. • In the ICM, each arc ( i, j ) ∈ A is assigned an activation probability π ij . Thepropagation process is started with a given seed set S . If node i has been activatedat the beginning of step t , then during step t , it will attempt to activate each ofits inactive neighbor node j with an independent probability π ij . Besides, for thosenodes that are successfully activated at step t , they will remain active and attemptto activate their inactive neighbor nodes in step t + 1. When no more inactive nodescan be activated, the diffusion process is terminated. • In the LTM, each arc ( i, j ) ∈ A is associated with a predefined weight b ij ≥ (cid:80) i : ( i,j ) ∈A b ij ≤ j ∈ V ; and each node i ∈ V selects a thresholdvalue θ i randomly chosen from [0 ,
1] before the propagation process. Let S t denotethe activated set in the beginning of step t ( S := S is the seed set). Then an inactivenode j ∈ V can be activated during step t if and only if (cid:80) i ∈S t b ij ≥ θ j , i.e. the totalcontribution of its active neighbors’ influence weights exceeds its threshold value θ j .In analogy to the ICM, if node j is successfully activated in some step, it will remainactive in the following steps, and the entire propagation process stops until no morenodes can be activated.Although the influence spread of the ICM or LTM is a stochastic process, Kempeet al. [17] converted it into a (discrete) deterministic process. More specifically,let Ω be the set of all possible scenarios of influence spread. Each scenario ω ∈ Ωcorresponds to a subgraph G ω = ( V , A ω ) of G (called a live-arc graph) with theprobability p ω (satisfying (cid:80) ω ∈ Ω p ω = 1). Here arc ( i, j ) ∈ A ω indicates that inscenario ω , if node i is activated in the influence spread, then node j must be ac-tivated by it. Let σ ω ( S ) represent the number of total activated nodes in G ω , i.e., σ ω ( S ) = |{ i ∈ V : ∃ a directed path from a node in S to i in G ω }| . Then the influencefunction µ ( S ) can equivalently be calculated by µ ( S ) = (cid:88) ω ∈ Ω p ω σ ω ( S ) . (2)Kempe et al. [17] showed that the IMPs under both ICM and LTM can be representedvia live-arc graphs. Indeed, under the ICM, each arc ( i, j ) is independently determinedto be live with probability π ij to construct a live-arc graph G ω . Hence, the probability of G ω is p ω = (cid:81) ( i,j ) ∈A ω π ij (cid:81) ( i,j ) ∈A\A ω (1 − π ij ). Under the LTM, for each node j ∈ V , itselects at most one incoming arc, for example ( i, j ), to be live with probability b ij , anddoes not select any arc with the probability 1 − (cid:80) i : ( i,j ) ∈A b ij . As a result, each live-arcgraph G ω under the LTM has a probability p ω = (cid:81) j ∈V I j , where I j = b ij if there exists i ∈ V such that ( i, j ) ∈ A ω ; and I j = 1 − (cid:80) ( i,j ) ∈A ω b ij , otherwise.4iven a social network G and a set of scenarios Ω, we next review the SMCLPformulation for the IMP [11]. First, for each scenario ω ∈ Ω and node i ∈ V , we denote R ( G ω , i ) as the reachability set of nodes that can activate node i in live-arc graph G ω (i.e., R ( G ω , i ) = { j ∈ V : ∃ a directed path from node j to node i in G ω } ). Then, foreach i ∈ V and ω ∈ Ω, we introduce binary variables y i and z ωi to denote whether node i is selected as a seed node and whether node i can be activated in scenario ω , respectively,i.e., y i = (cid:26) , if node i ∈ V is chosen as a seed node;0 , otherwise; z ωi = (cid:26) , if node i ∈ V can be activated in scenario ω ∈ Ω;0 , otherwise . Using the above notations, G¨uney et al. in [12] formulated the IMP as an SMCLP:max y , z (cid:88) ω ∈ Ω p ω (cid:88) i ∈V z ωi (3a)s.t. (cid:88) j ∈R ( G ω ,i ) y j ≥ z ωi , ∀ ω ∈ Ω , ∀ i ∈ V , (3b) (cid:88) i ∈V y i ≤ K, (3c) y i ∈ { , } , ∀ i ∈ V , (3d) z ωi ∈ { , } , ∀ ω ∈ Ω , ∀ i ∈ V . (3e)In this formulation, the objective function (3a) maximizes the expected number ofinfluenced nodes in the social network G . Constraints (3b) indicate that node i inscenario ω can be activated if and only if at least one of the nodes in its reachability set R ( G ω , i ) is chosen as a seed node. Constraint (3c) limits the cardinality of the set of seednodes up to K . Finally, constraints (3d) and (3e) restrict variables y and z to be binary.Unfortunately, formulation (3) is computational intractable for the network withrealistic dimensions due to the large number of scenarios (for the IMP under the ICMand LTM, the number of scenarios are all exponential). Hence, the Monte-Carlo samplingapproach is often used to approximate the influence diffusion process in which a reasonablesize of equiprobable scenarios set Ω (cid:48) ⊆ Ω is generated and the objective function in(3a) is replaced by (cid:80) ω ∈ Ω (cid:48) p (cid:48) ω (cid:80) i ∈V z ωi where p (cid:48) ω = 1 / | Ω (cid:48) | and hence (cid:80) ω ∈ Ω (cid:48) p (cid:48) ω = 1. Forsimplicity, we continue to use Ω and p ω to represent the set of sampling scenarios andthe probability of occurrence of scenario ω , respectively.However, the problem size of formulation (3) is still too large. Indeed, both numbersof variables and constraints are O ( | Ω ||V| ), making it difficult to solve the problem bystandard MIP solvers or the BD approach in [12], especially when the number of scenarios | Ω | or the sizes of graphs G ω is large. In the next section, we shall resolve this difficultyby proposing two new presolving methods to reduce the problem size. In this section, by exploiting the problem structure of formulation (3), we propose twopresolving methods and provide some theoretical analysis to show their strengths toreduce the problem size of formulation (3). Specifically, subsection 3.1 studies the SCNAwhich aggregates the nodes in each SCC, in a given live-arc graph, into a single node,and subsection 3.2 studies the INA which extends the idea of the SCNA into applyingisomorphic nodes aggregations among different live-arc graphs.5 .1 Strongly connected nodes aggregation
Figure 1:
The example live-arc graph which includes two SCCs, { } and { } . In this subsection, we present a presolving method based on a given live-arc graph. Tobegin with, we consider the example live-arc graph (corresponds to some scenario ω ) inFigure 1. In this graph, we note that there exists a directed path (arc) from node 1 tonode 2, and as a result, if node 1 is activated by some seed node, node 2 can also beactivated. Conversely, the fact that there exists a direct path (arc) from node 2 to node 1implies that if node 2 is activated by some seed node, node 1 can also be activated. Thismeans that either (i) nodes 1 and 2 can be simultaneously activated by some seed node;or (ii) neither of them can be activated. Consequently, z ω = z ω must hold in formulation(3). This reveals some redundancy in formulation (3) as it uses two variables z ω and z ω and two constraints in (3b) without considering z ω = z ω . Indeed, to simplify theproblem formulation, we can just remove variable z ω (or z ω ) and its associated constraintin (3b) and add the objective coefficient of variable z ω (or z ω ) into that of variable z ω (or z ω ). In general, we have the following statement. Proposition 3.1.
Given a live-arc graph G ω and its two strongly connected nodes i and i (i.e., there is a path from node i to node i in G ω and vice versa), we can set z ωi = z ωi in formulation (3) , and it does not change the optimal value. According to Proposition 3.1, we can remove variable z ωi (or z ωi ) and the associatedconstraint in (3b), and add the objective coefficient of variable z ωi (or z ωi ) into that ofvariable z ωi (or z ωi ) in the objective function. Notice that for a given SCC in a live-arcgraph, as all of its nodes are strongly connected, we can recursively apply the presolvingmethod in Proposition 3.1 until there remains only a single variable and a single constraintin (3b) associated with this SCC in formulation (3). The above presolving method iscalled SCNA in the following.1 2 : { , } { , , } Figure 2:
The compact live arc graph obtained by aggregating each SCC into asingle node in Figure 1. Notice that there are only two nodes 1 2 (corresponds toSCC { , } ) and 3 4 5 (corresponds to SCC { , , } ) and one arc between thesetwo nodes. After applying the SCNA, the IMP can be equivalently constructed based on the(weight) compact live-arc graphs. To be more specific, by aggregating each SCC of G ω into a single node with the weight being the size of the SCC, we can get a compact (anddirected acyclic) graph, denoted as ¯ G ω = ( ¯ V ω , ¯ A ω ) where each node v in ¯ V ω represents adistinct SCC in the original live-arc graph and each arc ( i, j ) ∈ ¯ A ω denotes that thereexists a path between the two SCCs corresponding to nodes i and j (see Figure 2 for an6xample of this transformation of the graph in Figure 1). For notational purpose, let SC ωv denote the SCC in G ω which is associated with node v ∈ ¯ V ω . As the reachabilitysets of the nodes inside a given SCC SC ωv are the same, we use the notation R ( G ω , SC ωv )to represent the reachability set of this SCC, which is equal to each R ( G ω , j ), j ∈ SC ωv .It follows immediately that R ( G ω , SC ωv ) = (cid:91) k ∈R ( ¯ G ω ,v ) SC ωk , (4)where R ( ¯ G ω , v ) is the reachability set of node v in the compact live-arc graph ¯ G ω . Thenconstraints (3b) reduce to (cid:88) k ∈R ( ¯ G ω ,v ) (cid:88) j ∈SC ωk y j ≥ z ωv , ∀ ω ∈ Ω , ∀ v ∈ ¯ V ω . Based on the above notations, the reduced SMCLP formulation can be written asmax y , z (cid:88) ω ∈ Ω p ω (cid:88) v ∈ ¯ V ω |SC ωv | z ωv (5a)s.t. (cid:88) k ∈R ( ¯ G ω ,v ) (cid:88) j ∈SC ωk y j ≥ z ωv , ∀ ω ∈ Ω , ∀ v ∈ ¯ V ω , (5b) (cid:88) i ∈V y i ≤ K, (5c) y i ∈ { , } , ∀ i ∈ V , (5d) z ωv ∈ { , } , ∀ ω ∈ Ω , ∀ v ∈ ¯ V ω . (5e) Theorem 3.2.
The numbers of variables z ωv and constraints (5b) in formulations (5) are equal to (cid:80) ω ∈ Ω | ¯ V ω | . Theorem 3.2 shows that the numbers of variables z ωv and associated constraints (5b)are equal to the number of the SCCs in the original live-arc graphs. This is much smallerthan those in formulation (3), which are equal to the number of nodes as in formulation(3). As a result, formulation (5) should be much more efficiently solvable than formulation(3), especially for the case that the numbers of SCCs are much smaller than the numbersof nodes in the live-arc graphs. In addition, the fact that formulation (5) is built onthe compact and directed acyclic live-arc graphs can significantly reduce the runtime ofgraph searching, which plays an important role in improving the performance of the BDalgorithm (see section 4 further ahead). The presolving method SCNA performs reductions on two nodes i and i in the samelive-arc graph G ω where nodes i and i are strongly connected, or equivalently, thereachability sets of nodes i and i are the same, i.e., R ( G ω , i ) = R ( G ω , i ). In thissubsection, we concentrate on extending the result into the isomorphic nodes amongdifferent live-arc graphs. As it has been previously mentioned, two nodes i and i intwo different live-arc graphs G ω and G η are called isomorphic if their reachability setsare the same, i.e., R ( G ω , i ) = R ( G η , i ) . (6)To begin with, we give the following observation stating that the values of variables z are determined by the values of variables y in formulation (3).7 bservation 3.3. There must exist an optimal solution ( ¯ y , ¯ z ) of formulation (3) suchthat ¯ z ωi = min (cid:88) j ∈R ( G ω ,i ) ¯ y j , , ∀ ω ∈ Ω , ∀ i ∈ V . (7)In particular, if node i ∈ V has no incoming arc in G ω in some ω ∈ Ω, i.e., R ( G ω , i ) = { i } is a singleton, the associated constraint in (3b) reduces to y i ≥ z ωi . By Observation 3.3,we can get z ωi = min { , y i } = y i . As a result, we can perform a reduction on formulation(3) by aggregating z ωi := y i and removing the associated constraint in (3b). We call thisreduction singleton node aggregation (SNA). Indeed, this is exactly the “P2” presolvingmethod proposed in G¨uney et al. [12].Using Observation 3.3, we can further derive reductions based on two isomorphicnodes in different live-arc graphs. Proposition 3.4.
Suppose that (6) holds for some scenarios ω, η ∈ Ω and nodes i , i ∈ V . Then setting z ωi = z ηi in formulation (3) does not change the optimal value.Proof. By Observation 3.3, there must exist an optimal solution ( ¯ y , ¯ z ) such that¯ z ωi = min , (cid:88) j ∈R ( G ω ,i ) ¯ y j and ¯ z ηi = min , (cid:88) j ∈R ( G η ,i ) ¯ y j . As nodes i and i are isomorphic, it follows from (6) that ¯ z ωi = ¯ z ηi . The proof iscomplete.From Proposition 3.4, variables z ωi and z ηi can be aggregated together in formulation(3). In other words, we may remove variable z ηi (or z ωi ) together with its associatedconstraint (cid:88) j ∈R ( G η ,i ) y j ≥ z ηi or (cid:88) j ∈R ( G ω ,i ) y j ≥ z ωi from formulation (3), and add the objective coefficient of variable z ηi (or z ωi ) into thatof variable z ωi ( z ηi ). The above presolving method is called INA in the following. Remark 3.5.
The SCNA can be regarded as a special case of the INA, which is restrictedto aggregating the isomorphic (strongly connected) nodes inside each live-arc graph.However, implementing the SCNA is equivalent to identifying all SCCs in all live-arcgraphs with a complexity of O ( (cid:80) ω ∈ Ω ( |V| + |A ω | )) , which is much faster than that ofimplementing the INA. The latter requires checking whether or not condition (6) holdsfor all four tuples ( ω, η, i , i ) with a complexity of O ( | Ω | |V| ) . It is interesting to ask that after the exhaustive applications of the presolving strategiesSNA and INA, how many variables z (or constraints (3b)) remain in the reduced problem.In the following, we answer this question by considering a special case of the IMP, whichis built on a one-way bipartite social network and under the LTM. We refer to [8, 26] forthe applications of one-way bipartite social networks and [35, 36] for some theoreticalresults of a variant of bipartite social network problem.More formally, let G B = ( M ∪ N , A ) be a given one-way bipartite network, where allarcs in A are from nodes in set M to those in set N . For each scenario ω ∈ Ω, denoteits live-arc graph as G ω B = ( M ∪ N , A ω ). It is important to note that under the LTM, ineach live-arc graph G ω B , each node i , i ∈ M , does not have any incoming arc and eachnode j , j ∈ N , has at most one incoming arc. Therefore, under the LTM, the reachabilityset of a node i ∈ M is R ( G ω B , i ) = { i } , and the reachability set of a node j ∈ N is R ( G ω B , j ) = (cid:26) { j } , if j has no incoming arc; { i, j } , if there exists some i ∈ M such that ( i, j ) ∈ A ω . (8)8ext, for each node i ∈ M ∪ N , we define a scenarios setΩ( i ) = { ω ∈ Ω : R ( G ω B , i ) = { i } , i.e., node i has no incoming arc in graph G ω B } , (9)and for each arc ( i, j ) ∈ A , we define another scenarios setΩ( i, j ) = { ω ∈ Ω : R ( G ω B , j ) = { i, j } , i.e., arc ( i, j ) is in graph G ω B } . (10)By definition, it follows that (cid:40) Ω( i ) = Ω for i ∈ M ;( ∪ i ∈M Ω( i, j )) ∪ Ω( j ) = Ω for j ∈ N . (11)Furthermore, we have the following:(i) by applying the SNA for each i ∈ M ∪ N , we can aggregate z ωi := y i for all ω ∈ Ω( i ),and remove the associated constraints in (3b);(ii) by applying the INA for each ( i, j ) ∈ A , we can aggregate all variables z ωj , ω ∈ Ω( i, j ),into a single variable, denoted as z ij , and remove the redundant constraints in (3b).Therefore, the reduced formulation is given bymax y , z (cid:88) i ∈M∪N s i y i + (cid:88) ( i,j ) ∈A c ij z ij (12a)s.t. y i + y j ≥ z ij , ∀ ( i, j ) ∈ A , (12b) (cid:88) i ∈M∪N y i ≤ K, (12c) y i ∈ { , } , ∀ i ∈ M ∪ N , (12d) z ij ∈ { , } , ∀ ( i, j ) ∈ A , (12e)where s i = (cid:80) ω ∈ Ω( i ) p ω for i ∈ M ∪ N and c ij = (cid:80) ω ∈ Ω( i,j ) p ω for ( i, j ) ∈ A . Obvi-ously, both numbers of constraints and variables are linear with the size of network G B (independent of the number of scenarios | Ω | ). Theorem 3.6.
Suppose that the IMP is built on a given scenarios of the one-way bipartitenetwork G B and under the LTM. Applying the presolving methods SNA and INA onformulation (3) , the numbers of variables and constraints in the reduced formulation (12) are equal to |M| + |N | + |A| and |A| + 1 , respectively. We next provide more analysis results for formulation (12). To proceed, we note thatusing (11) and (cid:80) ω ∈ Ω p ω = 1, we have following properties on the objective coefficientsof formulation (12). Remark 3.7. (i) s i = 1 for all i ∈ M ; and (ii) s j + (cid:80) i : ( i,j ) ∈A c ij = 1 for all j ∈ N . Theorem 3.8.
The LP relaxation of formulation (12) is tight. Moreover, formulation (12) can be solved in strongly polynomial-time.Proof.
Clearly, if K ≥ |M| + |N | , point ( y , z ) = ( e , e ) is optimal for formulation (12) andits LP relaxation, where e is an all ones vector with appropriate dimension. As a result,the statement follows. Therefore, in the following, we consider the case K < |M| + |N | .Let (¯ y , ¯ z ) be an optimal solution of the LP relaxation of formulation (12). If thereexists some i ∈ M and j ∈ N such that ¯ y i < y j >
0, then we can construct anew point (ˆ y , ˆ z ) as follows:ˆ y i = ¯ y i + ε, ˆ y j = ¯ y j − ε, and ˆ y i = ¯ y i for i ∈ M ∪ N \{ i , j } ;ˆ z ij = max { ¯ z ij − ε, } for ( i, j ) ∈ A and ˆ z ij = ¯ z ij for ( i, j ) ∈ A with j (cid:54) = j ,9here ε > y , ˆ z ) is also feasiblefor the LP relaxation of formulation (12). Moreover, point (ˆ y , ˆ z ) must be optimal since (cid:88) i ∈M∪N s i ˆ y i + (cid:88) ( i,j ) ∈A c ij ˆ z ij − (cid:88) i ∈M∪N s i ¯ y i + (cid:88) ( i,j ) ∈A c ij ¯ z ij = s i (¯ y i + ε ) + s j (¯ y j − ε ) + (cid:88) i : ( i,j ) ∈A c ij max { ¯ z ij − ε, }− s i ¯ y i − s j ¯ y j − (cid:88) i : ( i,j ) ∈A c ij ¯ z ij = s i ε − s j ε + (cid:88) i : ( i,j ) ∈A c ij max {− ε, − ¯ z ij }≥ s i ε − s j ε − (cid:88) i : ( i,j ) ∈A c ij ε = ε s i − s j − (cid:88) i : ( i,j ) ∈A c ij = 0 , where the last equality follows from Remark 3.7. Recursively applying the above argument,we will obtain an optimal solution (˜ y , ˜ z ) of the LP relaxation of formulation (12) fulfillingwith the following:( (cid:63) ) if ˜ y j > j ∈ N , then ˜ y i = 1 for all i ∈ M must be hold.Furthermore, we can, without loss of generality, assume the followings on point (˜ y , ˜ z ).1) (cid:80) i ∈M∪N ˜ y i = K . Otherwise, we can increase ˜ y i , with ˜ y i <
1, without decreasing theobjective value (as (cid:80) i ∈M∪N ˜ y i < K < |M| + |N | ).2) ˜ z ij = min { ˜ y i + ˜ y j , } for all ( i, j ) ∈ A . Otherwise, we can increase ˜ z ij , with ˜ z ij < min { ˜ y i + ˜ y j , } , without decreasing the objective value.Together with ( (cid:63) ) and 1), point (˜ y , ˜ z ) must satisfy the following(i) if K ≥ |M| , then ˜ y i = 1 for all i ∈ M ; and(ii) if K < |M| , then ˜ y j = 0 for all j ∈ N .We next prove the statement in the theorem by treating cases (i) and (ii) separately.(i) In this case, ˜ y i = 1 for all i ∈ M and by 2), ˜ z ij = 1 for all ( i, j ) ∈ A . Setting y i = 1 for all i ∈ M and z ij = 1 for all ( i, j ) ∈ A , the LP relaxation of formulation (12)reduces tomax y (cid:88) j ∈N s j y j + |M| + (cid:88) ( i,j ) ∈A c ij : (cid:88) j ∈N y j = K − |M| , y j ∈ [0 , , ∀ j ∈ N . (13)Then point { ˜ y j } j ∈N must be optimal to formulation (13). On the other hand, supposethat s j ≥ · · · ≥ s j |N| where { j , . . . , j |N | } = N . It is easy to show that point { y (cid:48) j } j ∈N is also optimal to formulation (13) where y (cid:48) j τ = 1 for τ = 1 , . . . , K − |M| and y (cid:48) j τ = 0otherwise. We next extend point { y (cid:48) j } j ∈N to a higher dimensional point ( y (cid:48) , z (cid:48) ) ∈{ , } |M| + |N | × { , } |A| by additionally setting y (cid:48) i = 1 for all i ∈ M and z (cid:48) ij = 1 for all( i, j ) ∈ A . Apparently, the 0-1 point ( y (cid:48) , z (cid:48) ) must also be optimal to formulation (12)and its LP relaxation. This implies that the LP relaxation of formulation (12) is tightand formulation (12) is strongly polynomial-time solvable for this case.(ii) In this case, ˜ y j = 0 for all j ∈ N and by 2), ˜ z ij = ˜ y i for all ( i, j ) ∈ A . Setting y i = 0 for all i ∈ N and ˜ z ij = ˜ y i for all ( i, j ) ∈ A , the LP relaxation of formulation (12)reduces tomax y (cid:88) i ∈M (cid:88) j : ( i,j ) ∈A c ij y i : (cid:88) i ∈M y i = K, y i ∈ [0 , , ∀ i ∈ M . (14)10hen point { ˜ y i } i ∈M must be optimal to formulation (13). Using the same argument incase (i), we can also prove the statement in this case.The tightness result in Theorem 3.6 sheds a useful insight on general (nonbipartite)networks. In particular, when the considered network contains lots of special structures(like one-way bipartite subgraphs), the LP relaxation of formulation (3) is likely tobe tight or nearly tight. Indeed, this is consistent to the empirical findings in G¨uney[11] in which the author observed that the gap between the optimal objective value offormulation (3) and that of its LP relaxation is very small for the IMP built on real-worldsocial networks. In this section, we shall discuss the implementation of the SCNA and INA. To proceedit, we first describe the algorithms for identifying the presolving reductions in subsection4.1 and then present the integration of them with the BD algorithm in subsection 4.2.
First, we note that identifying the reductions of the SCNA is easy since we only need toidentify all SCCs in each live-arc graph G ω = ( V , A ω ). This can be done in linear time O ( |V| + |A ω | ) using, for example, the Kosaraju’s algorithm [31]. After identifying allSCCs, we aggregate the nodes in each SCC into a (virtual) single node and transformlive-arc graph G ω into the compact and directed acyclic live-arc graph ¯ G ω = ( ¯ V ω , ¯ A ω ). Inthe following, we shall concentrate on the algorithm to detect the reductions of the INA.To identify the reductions of the INA, we need to first precompute and store thereachability sets of all nodes in all live-arc graphs G ω and then detect the pairs thatsatisfy condition (6). G¨uney et al. [12] computed the reachability set of each node i ineach (original) live-arc graph G ω , i.e., R ( G ω , i ), by applying a reverse breadth-first search(BFS) starting from node i . The associated computational complexity is O ( |V| + |A ω | ).Here we notice that it is much faster to compute R ( G ω , i ) based on the compact graph¯ G ω . Indeed, as it has been mentioned in subsection 3.1, the reachability sets of nodesinside a given SCC SC ωv are the same, where v is the corresponding node in the compactgraph ¯ G ω . Hence, to compute the reachability sets of the nodes in SCC SC ωv , we onlyneed to compute the reachability set R ( G ω , SC ωv ). To compute the latter one, we canapply a reverse BFS on the compact graph ¯ G ω and use the relation (4). Apparently,the associated complexity is O ( |V| + | ¯ A ω | ), which is much faster than the one thatapplies a reverse BFS on the (original) live-arc graph G ω , especially when | ¯ A ω | is muchsmaller than |A ω | . Furthermore, the fact that ¯ G ω is a directed acyclic graph (as itdoes not contain any SCC) allows to perform a topological ordering to further speedup the procedure of computing the reachability sets of all nodes in ¯ G ω . To be morespecific, topological ordering for the directed acyclic graph ¯ G ω is a linear ordering ofnodes such that for each arc ( v , v ) ∈ ¯ A ω , node v comes before v in the ordering.In our implementation, we traverse all nodes in ¯ G ω to compute their reachability setsaccording to the topological ordering. In other words, when computing the reachabilityset of node v , the reachability sets of nodes v (cid:48) , v (cid:48) ∈ N − ω ( v ), have been computed, where N − ω ( v ) = { v (cid:48) : ( v (cid:48) , v ) ∈ ¯ A ω } is the set of node v ’s neighborhood nodes. Therefore,to compute node v ’s reachability set, we only need to traverse its neighborhood nodesand use the relation R ( ¯ G ω , v ) = ∪ v (cid:48) ∈N − ω ( v ) R ( ¯ G ω , v (cid:48) ) ∪ { v } . This avoids performing awhole reverse BFS and hence generally accelerates the computation of R ( ¯ G ω , v ). Asfor the memory consumption, we note that it requires O ( |V| | Ω | ) memory to store thereachability sets of all nodes in all live-arc graphs G ω to implement the INA. Although11he SCNA can alleviate the memory consumption (as for each SCC, only a single node’sreachability set needs to be stored), it is still unrealistic to store all nodes’ reachabilitysets, especially for a large network or a high number of scenarios. To overcome thisdrawback, we restrict to storing those nodes’ reachability sets whose sizes are smallerthan or equal to a given parameter MaxReacSize . The rationale behind this strategyis that the nodes with smaller reachability sets are more likely to be isomorphic, asillustrated in our computational results (see subsection 5.4 further ahead).We next apply the INA by detecting the pairs whose reachability sets are the same,i.e., (cid:91) k ∈R ( ¯ G ω ,v ) SC ωk = (cid:91) k ∈R ( ¯ G ω ,v ) SC ω k (15)for some v ∈ ¯ V ω and v ∈ ¯ V ω . To do this, we follow [1] to use a hashing-based method.The basic idea of the hashing-based method is to simultaneously build a hashing tablethat remembers the information of the reachability sets and test whether there existsa reachability set in the hashing table that is identical to the one we are currentlylooking at. Specifically, let H be the hashing table and for each ω ∈ Ω and v ∈ ¯ V ω , let (cid:83) k ∈R ( ¯ G ω ,v ) SC ωk be the key with tuple ( ω, v, R ( ¯ G ω , v )) being the stored value in table H . At first, table H is initialized to be ∅ . Then, in each iteration, for scenario ω ∈ Ωand node v ∈ ¯ G ω (with | (cid:83) k ∈R ( ¯ G ω ,v ) SC ωk | ≤ MaxReacSize ), table H is queried for (cid:83) k ∈R ( ¯ G ω ,v ) SC ωk . If condition (15) holds for some corresponding entry ( ω , v , R ( ¯ G ω , v ))in table H , we apply the INA by removing variable z ωv , deleting the associated constraintin (5b), and adding the objective coefficient of variable z ωv into that of variable z ω v ;otherwise, tuple ( v, ω, R ( ¯ G ω , v )) will be added into table H . The procedure is repeateduntil all the considered reachability sets are tested.In summary, we present the implementation of the INA in Algorithm 1. Here,for simplicity of presentation, the improvement of topological ordering is omitted inAlgorithm 1. We use IsRemoved ( ω, v ) to indicate whether or not node v in scenario ω isaggregated by some node. In step 4, we perform the reverse BFS on node v in graph ¯ G ω to compute R ( ¯ G ω , v ). In steps 5-8, we use the hashing-based method to detect whetherthere exists some entry ( v , ω , R ( ¯ G ω , v )) in table H such that (15) holds and applythe INA reductions or add the new entry ( ω, v, R ( ¯ G ω , v )) into table H .After applying Algorithm 1, formulation (5) reduces tomax y , z (cid:88) ω ∈ Ω (cid:88) v ∈ ˜ V ω f ωv z ωv (16a)s.t. (cid:88) k ∈R ( ¯ G ω ,v ) (cid:88) j ∈SC ωk y j ≥ z ωv , ∀ ω ∈ Ω , ∀ v ∈ ˜ V ω , (16b) (cid:88) i ∈V y i ≤ K, (16c) y i ∈ { , } , ∀ i ∈ V , (16d) z ωv ∈ { , } , ∀ ω ∈ Ω , ∀ v ∈ ˜ V ω , (16e)where ˜ V ω = (cid:8) v ∈ ¯ V ω : IsRemoved ( ω, v ) = false (cid:9) . It deserves to mention that formu-lation (16) can be further simplified by applying the SNA: aggregating z ωv := y i andremoving the associated constraint in (3b) if (cid:83) k ∈R ( ¯ G ω ,v ) SC ωk = { i } . For simplicity ofpresentation, we omit this simplification in the following. G¨uney et al. [12] has proposed the BD algorithm to solve the IMP based on formulation(3). The author showed the effectiveness of integrating the presolving method SNA into12 lgorithm 1:
Implementation of the INA
Input: the compact live-arc graphs ¯ G ω = ( ¯ V ω , ¯ A ω ) and the SCCs SC ωv , v ∈ ¯ V ω ,in the (original) live-arc graphs G ω = ( V , A ω ), ω ∈ Ω. Output: the objective coefficients f ωv of variables z ωv and the indicatorsIsRemoved ( ω, v ), ω ∈ Ω and v ∈ ¯ V ω . Initialize IsRemoved ( ω, v ) = false , f ωv = p ω |SC ωv | ( ∀ ω ∈ Ω , v ∈ ¯ V ω ), and H = ∅ ; for ω ∈ Ω do for v ∈ ¯ V ω do Perform a reverse BFS on node v in graph ¯ G ω to compute R ( ¯ G ω , v ); if (cid:83) k ∈R ( ¯ G ω ,v ) |SC ω ( k ) | ≤ MaxReacSize and (15) holds for some ( v , ω , R ( ¯ G ω , v )) ∈ H then Set f ω v := f ω v + f ωv and IsRemoved ( ω, v ) := true ; else H := H (cid:83) { ( v, ω, R ( ¯ G ω , v )) } ; end end end the BD algorithm. In this subsection, to further enhance the capability of using the BDalgorithm to solve the IMP, we attempt to integrate the proposed SCNA and INA intothe BD algorithm, or equivalently, to design a BD algorithm that is based on the reducedformulation (16) and the compact graphs ¯ G ω , ω ∈ Ω.For notational convenience, for each ω ∈ Ω and v ∈ ¯ V ω , we denote (cid:80) j ∈SC ωv y j as I ωv ( y ), where SC ωv is the SCC (in the original live-arc graph G ω ) corresponding to node v (in the compact live-arc graph ¯ G ω ).4.2.1 Reformulation of (16)We first briefly introduce the reformulation of (16) based on the BD (more details canbe found in [12]). To begin with, we note that each binary variable z ωv can be relaxed tocontinuous variable since this does not change the optimal value of formulation (16). Let ϕ ω , ω ∈ Ω, represent the additional variable that captures the contribution of scenario ω to the objective function. Then, we can project out variables z and reformulate (16)equivalently as max ϕ , y (cid:88) ω ∈ Ω ϕ ω (17a)s.t. ϕ ω ≤ Φ ω ( y ) , ∀ ω ∈ Ω , (17b) (cid:88) i ∈V y i ≤ K, (17c) y i ∈ { , } , ∀ i ∈ V , (17d)where for ¯ y ∈ [0 , |V| and ω ∈ Ω, function Φ ω (¯ y ) is defined as follows:Φ ω (¯ y ) = max z (cid:88) v ∈ ˜ V ω f ωv z ωv : z ωv ≤ (cid:88) k ∈R ( ¯ G ω ,v ) I ωk (¯ y ) , ≤ z ωv ≤ , ∀ v ∈ ˜ V ω . (18)13or a fixed ¯ y ∈ [0 , |V| to model the inequality (17b), we use the Benders optimalitycuts which are derived as follows. First, the dual of formulation (18) ismin α ω , β ω (cid:88) v ∈ ˜ V ω α ωv (cid:88) k ∈R ( ¯ G ω ,v ) I ωk (¯ y ) + β ωv : α ωv + β ωv ≥ f ωv , α ωv , β ωv ≥ , ∀ v ∈ ˜ V ω , (19)where α ωv and β ωv are the dual variables of constraints z ωv ≤ (cid:80) k ∈R ( ¯ G ω ,v ) I ωk (¯ y ) and z ωv ≤ ¯ α ω (¯ y ) , ¯ β ω (¯ y )) where(¯ α ωv (¯ y ) , ¯ β ωv (¯ y )) = (cid:26) (0 , f ωv ) , if (cid:80) j ∈R ( ¯ G ω ,v ) I ωk (¯ y ) ≥ f ωv , , otherwise , ∀ v ∈ ˜ V ω . (20)Then the Benders optimality cuts for formulation (17) are given by ϕ ω ≤ (cid:88) v ∈ ˜ V ω ¯ α ωv (¯ y ) (cid:88) k ∈R ( ¯ G ω ,v ) I ωk ( y ) + ¯ β ωv (¯ y ) , ∀ ω ∈ Ω . (21)For simplicity, we shall abbreviate ¯ α ωv (¯ y ) and ¯ β ωv (¯ y ) as ¯ α ωv and ¯ β ωv , respectively. Let C ω = (cid:80) v ∈ ˜ V ω ¯ β ωv and c ωj = (cid:80) v ∈J ωj ¯ α ωv , j ∈ V , where J ωj = (cid:110) v ∈ ˜ V ω : j ∈ (cid:83) k ∈R ( ¯ G ω ,v ) SC ωk (cid:111) .Then Benders optimality cuts (21) can be alternatively written as ϕ ω ≤ (cid:88) j ∈V c ωj y j + C ω , ∀ ω ∈ Ω . (22)Given a point ¯ y ∈ [0 , |V| , it is interesting to ask whether or not applying the SCNAor INA changes the Benders optimality cuts (21), or equivalently, whether or not theBenders optimality cuts (21) based on formulations (3), (5), and (16) are equivalent. Proposition 4.1.
Given a point ¯ y ∈ [0 , |V| and a scenario ω ∈ Ω , (i) only applyingthe SCNA will not change the Benders optimality cut; (ii) applying the INA may changethe Benders optimality cut.Proof. The proof is relegated to the appendix.4.2.2 Solution approachWe now describe the BD algorithm to solve the Benders formulation (17) of the IMP.To implement the BD algorithm, we use a branch-and-Benders-cut approach in which abranch-and-bound tree is created and the Benders optimality cuts (17) are separated ateach node. Following [12, 34], we initially include the submodular cuts ϕ ω ≤ (cid:88) j ∈V Φ ω ( e j ) y j , ∀ ω ∈ Ω , (23)to formulation (17) in order to bound variables ϕ ω ( ω ∈ Ω), where e j ( j ∈ V ) is the unitvector with appropriate dimension.Having a solution (¯ y , ¯ ϕ ) of the LP relaxation of formulation (17), we next attempt toseparate the Benders optimality cuts (21). To do this, we need to compute all reachabilitysets of all nodes of all scenarios. However, as it has been mentioned in subsection 4.1, it isunrealistic to compute and store all reachability sets of all nodes of all scenarios in prioridue to the large memory consumption. Hence, similar to [12], we introduce parameter MaxLimPerScen to denote the maximally allowed memory consumption per scenario.14n particular, for each scenario, we store the reachability sets of nodes according totheir topological ordering in the compact live-arc graph until the memory consumptionreaches
MaxLimPerScen . As a result, when computing Benders optimality cut (21), if R ( ¯ G ω , v ) has been stored, we call it directly; otherwise, we perform a reverse BFS tocompute R ( ¯ G ω , v ) on the fly. Moreover, to further improve the efficiency of computingBenders optimality cut (21), we can omit to compute R ( ¯ G ω , v ) if (¯ α ωv , ¯ β ωv ) = (0 , f ωv ) isknown in priori. More specifically, let S ω = { k ∈ ¯ V ω : I ωk (¯ y ) ≥ } and R ( S ω ) be the setof nodes in ˜ V ω that can be reached from a node in S ω . Then, for each v ∈ R ( S ω ), wemust have (cid:80) k ∈R ( ¯ G ω ,v ) I ωk (¯ y ) ≥
1, and hence (¯ α ωv , ¯ β ωv ) = (0 , f ωv ). In summary, we presentthe separation of Benders optimality cuts in Algorithm 2. Algorithm 2:
Separation of Benders optimality cuts (21)
Input: the compact live-arc graphs ¯ G ω = ( ¯ V ω , ¯ A ω ), ω ∈ Ω, formulation (16), andpoint (¯ y , ¯ ϕ ). Output: the set C of Benders optimality cuts which are violated by point (¯ y , ¯ ϕ ). Initialize C = ∅ ; for ω ∈ Ω do Compute S ω := (cid:8) k ∈ ¯ V ω : I ωk (¯ y ) ≥ (cid:9) ; Perform a reverse BFS to compute R ( S ω ) := (cid:110) v ∈ ˜ V ω : v can be reached from a node in set S ω in graph ¯ G ω (cid:111) ; Initialize C ω := (cid:80) v ∈R ( S ω ) f ωv and c ωj := 0, j ∈ V ; for v ∈ ˜ V ω \ R ( S ω ) do If R ( ¯ G ω , v ) is not stored in the memory, perform a reverse BFS tocompute it ; if (cid:80) k ∈R ( ¯ G ω ,v ) I ωk (¯ y ) ≥ then Set C ω := C ω + f ωv ; else Set c ωj := c ωj + f ωv for all j ∈ (cid:83) k ∈R ( ¯ G ω ,v ) SC ωk ; end end if ¯ ϕ ω > (cid:80) j ∈V c ωj ¯ y j + C ω then Set C := C ∪ (cid:110) ϕ ω ≤ (cid:80) j ∈V c ωj y j + C ω (cid:111) ; end end In Algorithm 2, C ω and c ωj for each j ∈ V are used to keep track of the constant termand the coefficient of variable y j in Benders optimality cut (21). For each ω ∈ Ω, weinitialize C ω and c ωj in step 5 and then sequentially update them depending on whetheror not (cid:80) k ∈R ( ¯ G ω ,v ) I ωk (¯ y ) ≥ ϕ ω ≤ (cid:80) j ∈V c ωj y j + C ω is violated by point (¯ y , ¯ ϕ ), we add it into the set of violatedBenders optimality cuts C .It is worth emphasizing the computational efficiency of our separation algorithm forthe Benders optimality cuts over the one in [12]. First, in Algorithm 2, to computereachability sets R ( S ω ) and R ( ¯ G ω , v ), we perform reverse BFSes in the compact live-arcgraph ¯ G ω , which is much faster than that in G¨uney [12] where they performed the reverseBFS in the (original) live-arc graph G ω . Second, for each ω ∈ Ω, the number of variables z ωv is | ˜ V ω | which is much smaller than that in G¨uney [12], and as a result, fewer reverseBFSes will be performed in our separation algorithm.15 Computational results
In this section, we present the computational results to show the effectiveness of the pro-posed SCNA and INA. We use the BD algorithm in subsection 4.2, which is implementedin C++ linked with IBM ILOG CPLEX optimizer 12.7. The Benders cuts are addedusing CALLABLE LIBRARIES under the default settings of branch-and-cut frameworkof the CPLEX. The time limit is set to 14400 seconds, and all the experiments wereperformed on a cluster of Intel(R) Xeon(R) Gold 6140 CPU @ 2.30GHz computers. Onlya single core was used in our experiments.
Our benchmark data set consists of eight real-world social networks. Four of them havebeen used in [12, 34] (MSG, GNU, HEP, and ENRON), and the other networks are fromthe SNAP database (DEEZER, FACEBOOK, EPINIONS, and TWITTER). The latterones are large-scale networks which are used to test the performance of the SCNA andINA in large-scale cases. For the undirected networks (GNU, HEP, FACEBOOK, andDEEZER), we convert them into directed networks by adding two directed arcs ( i, j ) and( j, i ) for each edge ( i, j ). Table 1 summarizes the basic information of these networkswhere ρ = |A| / |V| is used to reflect the connectivity property of networks. Table 1:
Eight real-world social networks.Network |V| |A| ρ DescriptionMSG 1899 59835 31.5 Messaging network of UC-Irvine [27]GNU 10879 79988 7.4 Gnutella peer-to-peer file sharing network [29]HEP 15233 117782 7.7 High energy physics paper citation network [3]ENRON 36692 367662 10.0 Email communication network from Enron [21]FACEBOOK 50515 1638612 32.4 Facebook page network in the category of artist [30]DEEZER 54573 996404 18.3 Deezer friendship network of users in Croatia [30]TWITTER 81306 1768149 21.7 Social network from Twitter [19]EPINIONS 131828 841372 6.4 Who-trust-whom social network of Epinions [22]
We follow [12, 34] to choose the parameters in formulation (3) of the IMP as follows.For the IMP under both the ICM and LTM, the cardinality restriction K and the numberof scenarios | Ω | in formulation (3) are selected in { , , , } and { , , } ,respectively. For the IMP under the ICM, each arc ( i, j ) ∈ A is assigned the sameactivation probability, i.e., π ij = p , where p is chosen in { . , . , . } . For the IMPunder the LTM, we set influence weight b ij = 1 / deg( G , j ), where deg( G , j ) representsthe in-degree of node j in network G . In addition, one needs to take the effect of multipleparallel arcs in G into account. To be more specific, assume that there are n ij arcs fromnode i to node j . Under the ICM, π ij is set to 1 − (1 − p ) n ij to represent the probabilitythat at least one of these multiple parallel arcs ( i, j ) are activated. Under the LTM, b ij is set to n ij / deg( G , j ) as the influence weight needs to be counted n ij times. For theIMP with each combination of the above parameters, 5 instances are randomly generated.Therefore, for each social network in Table 1, we have 180 and 60 instances for the IMPunder the ICM and LTM, respectively.In our experiments, we compare the performance of the following three settings: • Default : solving the IMP based on the Benders reformulation of (3) with the SNAapplied; • SCNA : Default with the SCNA applied;16
INA : Default with the SCNA and INA applied.In our implementation, following [12], the memory control parameter
MemLimPerScen is set to 8 / | Ω | GB. Unless otherwise stated, in the implementation of the INA, parameter
MaxReacSize is set to 8 and 4 for the IMP under the ICM and LTM, respectively.Finally, to avoid generating too many Benders optimality cuts at fractional points, theseparation procedure stops if the dual bound improves by less than 0.001.
In this subsection, we test the effectiveness of the proposed presolving methods SCNAand INA for the IMP under the ICM.We first compare the sizes of the compact live-arc graphs ¯ G ω = ( ¯ V ω , ¯ A ω ) (obtainedby applying the SCNA) and the original live-arc graphs G ω = ( V , A ω ). Table 2 lists theaverage percentages of nodes and arcs reductions in 1000 scenarios (∆V = ( |V|−| ¯ V ω | ) / |V| and ∆A = ( |A ω | − | ¯ A ω | ) / |A ω | ). From the table, we can clearly observe that with theSCNA, the numbers of nodes and arcs in the compact live-arc graphs are much smallerthan those in the original live-arc graphs. In addition, we can see that the larger theactivation probability p is, the more reductions that will be detected by the SCNA. Thisis reasonable since as p increases, the original live-arc graphs contain more arcs. As aresult, more and larger SCCs are likely to appear. The same argument is applied tothe case that the network has a large connectivity ρ . For instance, the SCNA is moreeffective for networks MSG and FACEBOOK as they have a relative large value ρ . Table 2:
Average reductions of nodes and arcs for the IMP under the ICMthrough applying the SCNA.Network p ∆V ∆A Network p ∆V ∆AMSG 0.01 2.2% 12.9% FACEBOOK 0.01 2.4% 13.9%( ρ = 31 .
5) 0.05 25.6% 78.8% ( ρ = 32 .
4) 0.05 31.7% 77.0%0.10 37.2% 88.6% 0.10 51.4% 90.7%GNU 0.01 0.0% 1.0% DEEZER 0.01 0.1% 1.0%( ρ = 7 .
4) 0.05 0.9% 5.1% ( ρ = 18 .
3) 0.05 16.6% 39.1%0.10 8.3% 19.9% 0.10 44.8% 77.1%HEP 0.01 0.2% 4.6% TWITTER 0.01 2.1% 24.4%( ρ = 7 .
7) 0.05 4.9% 29.7% ( ρ = 21 .
7) 0.05 18.7% 67.5%0.10 13.4% 52.1% 0.10 34.2% 83.3%ENRON 0.01 0.1% 2.2% EPINIONS 0.01 0.2% 4.1%( ρ = 10 .
0) 0.05 7.6% 50.0% ( ρ = 6 .
4) 0.05 3.6% 52.5%0.10 15.5% 66.1% 0.10 6.4% 65.9%
Next, we report the reductions of variables and nonzeros in formulation (3) byapplying the SNA, SCNA, and INA in Table 3. We use ∆VAR and ∆NNZ to representthe percentage of the eliminated variables and nonzeros in formulation (3) throughapplying the SNA (which has been embedded in setting
Default ). In addition, we use+∆VAR and +∆NNZ to denote the additional eliminated variables and nonzeros (beyondthe SNA) through applying the SCNA or INA. As it can be seen in Table 3, when thenetwork has a small connectivity ρ or activation probability p , the singleton nodes (nodeswithout incoming arc) are more likely to appear in the live-arc graphs. As a result, theSNA can eliminate a considerably large number of variables and nonzeros. In contrast,the SCNA is more effective in eliminating the variables and nonzeros when the network17as a relative large value of ρ or p since as it has been previously mentioned, more andlarger SCCs are more likely to appear in these cases. For the INA, since the isomorphicnodes inside each scenario has been identified by the SCNA, the additional reductionscomes from detecting the isomorphic nodes among different scenarios. Nevertheless, wecan observe a clear reduction on the number of variables (beyond the SCNA). However,the reduction on the number of nonzeros is relatively small, which is due to the factthat most eliminated variables are associated with small reachability sets as we set MaxReacSize = 8 (in subsection 5.4, we will perform numerical experiments to confirmthat setting
MaxReacSize = 8 is enough to identify most pairs of isomorphic nodes).
Table 3:
Average reductions of the variables and nonzeros for the IMP underthe ICM through applying the SNA, SCNA, and INA.
Default SCNA INA
Network p ∆VAR ∆NNZ +∆VAR +∆NNZ +∆VAR +∆NNZMSG 0.01 80.9% 16.1% +2.1% +19.4% +6.9% +20.9%( ρ = 31 .
5) 0.05 56.4% 0.4% +25.6% +62.4% +27.9% +62.4%0.10 44.1% 0.2% +37.2% +68.2% +38.6% +68.3%GNU 0.01 93.1% 89.3% +0.0% +0.1% +5.4% +7.9%( ρ = 7 .
4) 0.05 72.3% 46.2% +0.9% +2.3% +14.5% +16.2%0.10 55.8% 0.3% +8.3% +36.1% +19.0% +36.2%HEP 0.01 93.3% 88.4% +0.2% +0.4% +5.3% +8.0%( ρ = 7 .
7) 0.05 75.8% 3.1% +4.9% +49.3% +17.8% +50.2%0.10 62.1% 0.4% +13.4% +65.4% +29.7% +65.6%ENRON 0.01 92.8% 34.4% +0.1% +6.5% +3.1% +8.1%( ρ = 10 .
0) 0.05 76.8% 0.2% +7.6% +48.8% +13.6% +48.8%0.10 64.1% 0.1% +15.5% +56.6% +22.8% +56.6%FACEBOOK 0.01 79.1% 0.6% +2.4% +31.0% +8.5% +31.0%( ρ = 32 .
4) 0.05 47.8% 0.0% +31.7% +69.6% +36.3% +69.6%0.10 32.5% 0.0% +51.4% +80.4% +54.3% +80.4%DEEZER 0.01 84.4% 73.7% +0.1% +0.2% +9.2% +12.2%( ρ = 18 .
3) 0.05 50.6% 0.0% +16.6% +54.0% +26.5% +54.0%0.10 32.4% 0.0% +44.8% +75.8% +50.6% +75.8%TWITTER 0.01 83.8% 0.9% +2.1% +44.9% +8.0% +45.0%( ρ = 21 .
7) 0.05 57.9% 0.0% +18.7% +64.7% +26.2% +64.7%0.10 43.1% 0.0% +34.2% +72.4% +40.6% +72.4%EPINIONS 0.01 95.9% 19.4% +0.2% +16.1% +1.7% +16.5%( ρ = 6 .
4) 0.05 87.7% 0.2% +3.6% +44.5% +6.8% +44.5%0.10 81.2% 0.1% +6.4% +46.5% +10.8% +46.5%
We now evaluate the performance improvement of the integration of the SCNA andINA with the BD algorithm. For each of the eight networks, Table 4 reports the totalnumber of instances that can be solved within 14400 seconds (
S), the average CPUtime in seconds (T), the average presolving time in seconds (PT), the average numberof branch nodes (
N), and the average number of added Benders optimality cuts (
C).Notice that the average CPU time T includes the presolving time spent on applying theSCNA and INA. Detailed statistics of these results can be found in Tables 10a-17a inthe appendix.As it can be seen in Table 4, the performance of setting
SCNA is much better18han that of setting
Default , especially for instances with large-scale networks. Intotal, setting
SCNA can solve 1248 instances with an average CPU time 414 . Default can only solve 1070 instances with an average CPU time 821 . Default and
SCNA . Notice that it isreasonable to observe that the numbers of added cuts in other networks are differentsince some instances cannot be solved within the time limit.As for setting
INA , we observe that it slightly improves the performance, comparedto setting
SCNA . In total, setting
INA can solve 13 more instances than setting
SCNA ,with the average CPU time decreasing from 414 .
0s to 332 . Table 4:
Performance improvement for the IMP under the ICM through applyingthe SCNA and INA.
Default SCNA INA
Network
S T N C S T PT N C S T PT N CMSG
180 2.7
170 78.8
176 141.0
164 517.9
106 3427.8
120 1635.6
165 1587.7
180 778.2
In this subsection, we present similar computational results for the IMP under the LTMin Tables 5-7. To begin with, we note that for a node in the original graph G , if it is asingleton node, (i.e., it does not have any incoming arc), then in each live-arc graph G ω ,it is also a singleton node; otherwise, it has exactly one incoming arc in each live-arcgraph G ω (due to the fact that we choose b ij = 1 / deg( G , j ) for all ( i, j ) ∈ A in our testedinstances). As a result, there exists at most one path between every pair of nodes in eachlive-arc graph G ω , implying that each SCC must be a circle (which is in sharp contrastto the IMP under the ICM in which each SCC can be the union of multiple circles). Thisproperty for the IMP under the LTM, however, implies that in the live-arc graphs, theSCCs are likely to be small, and (hence) the reachability sets of the nodes are likely to besmall. Consequently, through applying the SCNA, we can only observe a mild reductionon the numbers of nodes and arcs in Table 5 as well as a slight reduction on the numbers19f variables and nonzeros in Table 6. However, due to the small sizes of the reachabilitysets, more nodes are likely to be isomorphic among different scenarios and hence the INAcan detected more reductions, as compared to the IMP under the ICM. This is shownin Table 6 in which we observe a large reduction on the number of variables (beyondapplying the SCNA). Similarly, the reduction on the number of nonzeros is relativelysmall which can be explained by the fact that most eliminated variables by the INA areassociated with small reachability sets as MaxReacSize is set to 4 in our implementation(see subsection 5.4 further ahead for the reason of setting
MaxReacSize = 4). Noticethat in Table 6, the reduction detected by the SNA (
Default ) is marginal since, as ithas been mentioned, a node is a singleton node in live-arc graph G ω if and only if it is asingleton node in the original graph G . Therefore, the reduction of the SNA for the IMPunder the LTM totally depends on the number of singleton nodes in network G , whichis very small in most cases. Indeed, only network EPINIONS contains a relative largepercentage of singleton nodes (35 . Table 5:
Average reductions of nodes and arcs for the IMP under the LTMthrough applying SCNA.Network ∆V ∆A Network ∆V ∆AMSG 2.4% 4.6% FACEBOOK 2.9% 5.5%GNU 5.9% 11.8% DEEZER 4.4% 8.5%HEP 22.7% 42.6% TWITTER 3.3% 5.4%ENRON 8.7% 16.0% EPINIONS 1.8% 5.5%
Table 6:
Average reductions of the variables and nonzeros for the IMP underthe LTM through applying the INA, SCNA, and INA.
Default SCNA INA
Network ∆VAR ∆NNZ +∆VAR +∆NNZ +∆VAR +∆NNZMSG 1.9% 0.3% +2.4% +0.8% +12.6% +4.5%GNU 0.0% 0.0% +5.9% +1.9% +26.2% +10.5%HEP 0.0% 0.0% +22.7% +15.3% +71.7% +54.2%ENRON 0.0% 0.0% +8.7% +3.0% +31.3% +12.1%FACEBOOK 0.0% 0.0% +2.9% +0.4% +7.2% +0.9%DEEZER 0.0% 0.0% +4.4% +0.9% +11.6% +2.4%TWITTER 0.0% 0.0% +3.3% +1.0% +8.4% +2.2%EPINIONS 35.9% 10.3% +1.8% +0.8% +21.8% +10.8%
We now present the overall performance improvement of integrating the SCNA andINA into the BD algorithm in Table 7. Detailed statistics of these results can be foundin Tables 10b-17b in the appendix. From Table 7, we can see that the performanceof settings
SCNA and
INA is slightly better than setting
Default . Indeed, we onlyobserve a minor reduction on the average runtime (T) and the number of solved instances(
S) through applying the SCNA and INA. This can be explained by the fact that (i) thereduction on the sizes of networks through applying SCNA is small (as shown in Table 5)and (ii) only isomorphic nodes with small reachability sets can be detected by the INA(as shown in column +∆NNZ of setting
INA in Table 6). In addition, we note fromTable 7 that the IMPs under the LTM are much easier than those under the ICM. Intotal, under the LTM, only 20 among 480 instances (4.2%) cannot be solved by setting
INA within the given time limit while 179 among 1440 instances (12.4%) cannot be20olved by the same setting under the ICM.
Table 7:
Performance improvement for the IMP under the LTM through applyingthe SCNA and INA.
Default SCNA INA
Network
S T N C S T PT N C S T PT N CMSG
60 23.7
60 62.8
60 55.6
60 203.6
47 1689.3
60 334.0
60 724.7
460 242.3
As it has been mentioned in subsection 4.1,
MaxReacSize is a parameter to achieve atrade-off between the effectiveness and efficiency of implementing the INA: the largerthe parameter
MaxReacSize is, the more isomorphic nodes that will be identified andthe higher computational complexity. Therefore, in this subsection, we compare theperformance of different selections of parameter
MaxReacSize . Tables 8 and 9 reportthe computational results for the IMP under the ICM and LTM, respectively (for theICM, we only report the results of the case p = 0 . p = 0 .
01 or 0 . to represent the total memory consumption (in GB) of storing all reachabilitysets after removing those detected by the SCNA. In addition, we use M to denote thememory consumption (in GB) of storing the reachability sets with the size restriction( MaxReacSize = 2 , , , , INA to denote the associated computational time ofimplementing INA, and ∆VAR
INA to denote the average extra percentage of variablesthat can be eliminated by the INA.For the ICM, Table 8 shows that it requires a prohibitively large memory M to storeall the reachability sets. However, when restricting the size of the considered reachabilitysets to a small value MaxReacSize , the memory overhead significantly reduces; seecolumn M in Table 8. In addition, with the increasing value of
MaxReacSize , theimprovement on the percentage of the eliminated variables ∆VAR
INA becomes smallerand smaller. Indeed, when
MaxReacSize ≥
8, ∆VAR
INA almost keeps unchanged.This shows that with a small value of parameter
MaxReacSize , our algorithm canidentify almost all pairs of isomorphic nodes. We note that for the IMP under theICM, the computational time of implementing the INA is relatively small, even when
MaxReacSize = 10.We now discuss the results for the IMP under the LTM in Table 9. On one hand, forthe IMP under the LTM, since there exists at most one incoming arc for each node ineach live-arc graph, the sizes of the reachability sets are likely to be smaller than thosefor the IMP under the ICM. This leads to a smaller total memory consumption M anda larger memory consumption M when restricting MaxReacSize to a small value, ascompared to those for the IMP under the ICM. As a result, the computational overheadof implementing the INA is very high, especially when parameter
MaxReacSize is large.In addition, in analogy to the IMP under the ICM, with a small value of parameter21 a b l e : C o m p a r i s o n o f e ff e c t i v e n e ss o f t h e I NA w i t hd i ff e r e n t s e l e c t i o n s o f M a x R e a c S i z e ( f o r t h e I M P und e r t h e I C M ) . M a x R e a c S i z e N e t w o r k M M T I NA ∆ VA R I NA M T I NA ∆ VA R I NA M T I NA ∆ VA R I NA M T I NA ∆ VA R I NA M T I NA ∆ VA R I NA M S G . . . . % . . . % . . . % . . . % . . . % G NU . . . . % . . . % . . . % . . . % . . . % H EP . . . . % . . . % . . . % . . . % . . . % E N R O N . . . . % . . . % . . . % . . . % . . . % F A C E B OO K . . . . % . . . % . . . % . . . % . . . % D EE Z E R . . . . % . . . % . . . % . . . % . . . % T W I TT E R . . . . % . . . % . . . % . . . % . . . % EP I N I O N S . . . . % . . . % . . . % . . . % . . . % T a b l e : C o m p a r i s o n o f e ff e c t i v e n e ss o f t h e I NA w i t hd i ff e r e n t s e l e c t i o n s o f M a x R e a c S i z e ( f o r t h e I M P und e r t h e L T M ) . M a x R e a c S i z e N e t w o r k M M T I NA ∆ VA R I NA M T I NA ∆ VA R I NA M T I NA ∆ VA R I NA M T I NA ∆ VA R I NA M T I NA ∆ VA R I NA M S G . . . . % . . . % . . . % . . . % . . . % G NU . . . . % . . . % . . . % . . . % . . . % H EP . . . . % . . . % . . . % . . . % . . . % E N R O N . . . . % . . . % . . . % . . . % . . . % F A C E B OO K . . . . % . . . % . . . % . . . % . . . % D EE Z E R . . . . % . . . % . . . % . . . % . . . % T W I TT E R . . . . % . . . % . . . % . . . % . . . % EP I N I O N S . . . . % . . . % . . . % . . . % . . . % axReacSize , our algorithm can identified almost all isomorphic nodes for the IMPunder the LTM. Therefore, for the IMP under the LTM, we choose MaxReacSize = 4in the implementation of the INA to achieve a trade-off between the performance andthe time complexity.
In this paper, we have proposed two new presolving methods, called the SCNA and INA,and integrated them into the BD algorithm to solve the IMP. The SCNA enables to buildan SMCLP formulation for the considered problem based on the compact live-arc graphs,which are obtained by aggregating strongly connected nodes in the original live-arc graphs.The INA further reduces the problem size of the SMCLP formulation. In particular, forthe IMP built on the one-way bipartite network G B under the LTM, through applyingthe INA, we are able to (i) provide a compact formulation whose problem size dependsonly on the size of network G B but not on the number of the Monte-Carlo sampling;and (ii) show that this compact formulation is strongly polynomial-time solvable and itsLP relaxation is tight. Furthermore, with the SCNA and INA, we can develop a muchfaster separation procedure for the Benders optimality cuts, which plays a crucial rolein speeding up the BD algorithm. We perform extensive experiments to analyze theperformance impact of the proposed SCNA and INA on solving the IMP with real-worldsocial networks. Computational results show that the proposed SCNA and INA cansignificantly reduce the problem size and hence improve the performance of using the BDalgorithm to solve the IMP, especially for problems with large-scale networks and underthe ICM. For future work, it is interesting to develop more powerful presolving methodsfor solving the IMP. In addition, it also deserves to investigate whether the proposedpresolving methods are effective in solving other variants of IMPs [35, 36]. References [1] T. Achterberg, R. E. Bixby, Z. Gu, E. Rothberg, and D. Weninger. Presolve reductions inmixed integer programming.
INFORMS Journal on Computing , 32(2):473–506, 2020.[2] R. Bornd¨orfer.
Aspects of set packing, partitioning and covering . PhD thesis, TechnischeUniversit¨at Berlin, 1998.[3] W. Chen, Y. Wang, and S. Yang. Efficient influence maximization in social networks. In
Proceedings of the 15th ACM SIGKDD International Conference on Knowledge Discoveryand Data Mining , pages 199–208, 2009.[4] W. Chen, C. Wang, and Y. Wang. Scalable influence maximization for prevalent viralmarketing in large-scale social networks. In
Proceedings of the 16th ACM SIGKDDInternational Conference on Knowledge Discovery and Data Mining , pages 1029–1038,2010.[5] W. Chen, Y. Yuan, and L. Zhang. Scalable influence maximization in social networks underthe linear threshold model. In , pages88–97, 2010.[6] S. Cheng, H. Shen, J. Huang, G. Zhang, and X. Cheng. Staticgreedy: Solving the scalability-accuracy dilemma in influence maximization. In
Proceedings of the 22nd ACM InternationalConference on Information & Knowledge Management , pages 509–518, 2013.[7] CPLEX. . 2020.[8] G. Erg¨un. Human sexual contact network as a bipartite graph.
Physica A: StatisticalMechanics and its Applications , 308(1):483–488, 2002.
9] M. Fischetti, M. Kahr, M. Leitner, M. Monaci, and M. Ruthmair. Least cost influencepropagation in (social) networks.
Mathematical Programming , 170(1):293–325, 2018.[10] S. Galhotra, A. Arora, and S. Roy. Holistic influence maximization: Combining scalabil-ity and efficiency with opinion-aware models. In
Proceedings of the 2016 InternationalConference on Management of Data , pages 743–758, 2016.[11] E. G¨uney. An efficient linear programming based method for the influence maximizationproblem in social networks.
Information Sciences , 503:589–605, 2019.[12] E. G¨uney, M. Leitner, M. Ruthmair, and M. Sinnl. Large-scale influence maximization viamaximal covering location.
European Journal of Operational Research , 289(1):144 – 164,2021.[13] D. G¨unne¸c, S. Raghavan, and R. Zhang. Least-cost influence maximization on socialnetworks.
INFORMS Journal on Computing , 32(2):289–302, 2020.[14] GUROBI. . 2020.[15] S. Heinz, J. Schulz, and J. C. Beck. Using dual presolving reductions to reformulatecumulative constraints.
Constraints , 18(2):166–201, 2013.[16] M. Kahr, M. Leitner, M. Ruthmair, and M. Sinnl. Benders decomposition for competitiveinfluence maximization in (social) networks.
Omega, accepted for publication , 2020.[17] D. Kempe, J. Kleinberg, and ´E. Tardos. Maximizing the spread of influence through asocial network. In
Proceedings of the 9th ACM SIGKDD International Conference onKnowledge Discovery and Data Mining , pages 137–146, 2003.[18] M. Kimura and K. Saito. Tractable models for information diffusion in social networks. In
Knowledge Discovery in Databases: PKDD 2006 , pages 259–271, 2006.[19] J. Leskovec and J. J. Mcauley. Learning to discover social circles in ego networks. In
Advances in Neural Information Processing Systems , pages 539–547, 2012.[20] J. Leskovec, A. Krause, C. Guestrin, C. Faloutsos, J. VanBriesen, and N. Glance. Cost-effective outbreak detection in networks. In
Proceedings of the 13th ACM SIGKDDInternational Conference on Knowledge Discovery and Data Mining , pages 420–429, 2007.[21] J. Leskovec, K. J. Lang, A. Dasgupta, and M. W. Mahoney. Community structure in largenetworks: Natural cluster sizes and the absence of large well-defined clusters.
InternetMathematics , 6(1):29–123, 2009.[22] J. Leskovec, D. Huttenlocher, and J. Kleinberg. Signed networks in social media. In
Proceedings of the SIGCHI Conference on Human Factors in Computing Systems , pages1361–1370, 2010.[23] X. Li, J. D. Smith, T. N. Dinh, and M. T. Thai. Tiptop: (almost) exact solutions forinfluence maximization in billion-scale networks.
IEEE/ACM Transactions on Networking ,27(2):649–661, 2019.[24] Y. Li, J. Fan, Y. Wang, and K.-L. Tan. Influence maximization on social graphs: A survey.
IEEE Transactions on Knowledge and Data Engineering , 30(10):1852–1872, 2018.[25] I. Ljubi´c, P. Putz, and J.-J. Salazar-Gonz´alez. Exact approaches to the single-sourcenetwork loading problem.
Networks , 59(1):89–106, 2012.[26] M. E. J. Newman. The structure and function of complex networks.
SIAM Review , 45(2):167–256, 2003.[27] P. Panzarasa, T. Opsahl, and K. M. Carley. Patterns and dynamics of users’ behavior andinteraction: Network analysis of an online community.
Journal of the American Society forInformation Science and Technology , 60(5):911–932, 2009.
28] S. Raghavan and R. Zhang. A branch-and-cut approach for the weighted target set selectionproblem on social networks.
INFORMS Journal on Optimization , 1(4):304–322, 2019.[29] M. Ripeanu and I. Foster. Mapping the gnutella network: Macroscopic properties oflarge-scale peer-to-peer systems. In
International Workshop on Peer-to-Peer Systems , pages85–93, 2002.[30] B. Rozemberczki, R. Davies, R. Sarkar, and C. Sutton. Gemsec: Graph embedding with selfclustering. In
Proceedings of the 2019 IEEE/ACM International Conference on Advancesin Social Networks Analysis and Mining 2019 , pages 65–72, 2019.[31] M. Sharir. A strong-connectivity algorithm and its applications in data flow analysis.
Computers & Mathematics with Applications , 7(1):67–72, 1981.[32] Y. Tang, X. Xiao, and Y. Shi. Influence maximization: Near-optimal time complexity meetspractical efficiency. In
Proceedings of the 2014 ACM SIGMOD International Conference onManagement of Data , pages 75–86, 2014.[33] Y. Tang, Y. Shi, and X. Xiao. Influence maximization in near-linear time: A martin-gale approach. In
Proceedings of the 2015 ACM SIGMOD International Conference onManagement of Data , page 1539–1554, 2015.[34] H.-H. Wu and S. K¨u¸c¨ukyavuz. A two-stage stochastic programming approach for influencemaximization in social networks.
Computational Optimization and Applications , 69(3):563–595, 2018.[35] H.-H. Wu and S. K¨u¸c¨ukyavuz. Probabilistic partial set covering with an oracle for chanceconstraints.
SIAM Journal on Optimization , 29(1):690–718, 2019.[36] P. Zhang, W. Chen, X. Sun, Y. Wang, and J. Zhang. Minimizing seed set selectionwith probabilistic coverage guarantee in a social network. In
Proceedings of the 20thACM SIGKDD International Conference on Knowledge Discovery and Data Mining , pages1306–1315, 2014.
A AppendixProof of Proposition 4.1
Proof.
Note that if only the SCNA is applied, formulation (16) reduces to formulation (5);and if neither the SCNA nor the INA is applied, formulation (16) reduces to formulation (3).Hence, we shall prove case (i) by showing that the Benders optimality cut based on formulation(3) is equivalent to the one based on formulation (5). First, for formulation (3), the Bendersoptimality cut is given by ϕ ω ≤ (cid:88) i ∈V ¯ α ωi (cid:88) j ∈R ( G ω ,i ) y j + ¯ β ωi (24)where for each i ∈ V , (¯ α ωi , ¯ β ωi ) = (0 , p ω ) if (cid:80) j ∈R ( G ω ,i ) ¯ y j ≥
1; and (¯ α ωi , ¯ β ωi ) = ( p ω ,
0) otherwise.Note that for all nodes in a same SCC SC ωv , their reachability sets are the same and equal to R ( G ω , SC ωv ). As a result, (¯ α ωi , ¯ β ωi ) for all i ∈ SC ωv must also be the same, denoted by (¯ α ωv , ¯ β ωv ).Hence, inequality (24) can be rewritten as ϕ ω ≤ (cid:88) v ∈ ¯ V ω |SC ωv | ¯ α ωv (cid:88) j ∈R ( G ω , SC ωv ) y j + |SC ωv | ¯ β ωv (25)As for formulation (5), the Benders optimality cut reads ϕ ω ≤ (cid:88) v ∈ ¯ V ω ¯ λ ωv (cid:88) k ∈R ( ¯ G ω ,v ) (cid:88) j ∈SC ωk y j + ¯ η ωv , (26) here for each v ∈ ¯ V , (¯ λ ωv , ¯ η ωv ) = (0 , p ω |SC ωv | ) if (cid:80) k ∈R ( ¯ G ω ,v ) (cid:80) j ∈SC ωk y j ≥
1; and (¯ λ ωv , ¯ η ωv ) =( p ω |SC ωv | ,
0) otherwise. Together with equation (4), we have (¯ λ ωv , ¯ η ωv ) = ( |SC ωv | ¯ α ωv , |SC ωv | ¯ β ωv ), andhence the two Benders optimality cuts (25) and (26) are equivalent. This proves the case (i).We next prove case (ii). Indeed, for two scenarios ω and η , the corresponding Bendersoptimality cuts based on formulation (5) are given by ϕ ω ≤ (cid:88) v ∈ ¯ V ω ¯ α ωv (cid:88) k ∈R ( ¯ G ω ,v ) I ωk ( y ) + ¯ β ωv , (27) ϕ η ≤ (cid:88) v ∈ ¯ V η ¯ α ηv (cid:88) k ∈R ( ¯ G η ,v ) I ηk ( y ) + ¯ β ηv . (28)Suppose that node v in graph ¯ G ω is isomorphic to node v in graph ¯ G η . By applying theINA, we can remove the variable z ωv and the associated constraint (5b) in formulation (5).The new objective coefficient f ηv is set to the sum of the old objective coefficients of variables z ωv and z ηv , while the new objective coefficient f ωv is set to zero. As a result, the term¯ α ωv (cid:80) k ∈R ( ¯ G ω ,v ) I ωk ( y ) + ¯ β ωv in Benders optimality cut (27) will be incorporated into Bendersoptimality cut (28), leading to two new Benders optimality cuts ϕ ω ≤ (cid:88) v ∈ ¯ V ω \{ v } ¯ α ωv (cid:88) k ∈R ( ¯ G ω ,v ) I ωk ( y ) + ¯ β ωv , (29) ϕ η ≤ (cid:88) v ∈ ¯ V η ¯ α ηv (cid:88) k ∈R ( ¯ G η ,v ) I ηk ( y ) + ¯ β ηv + ¯ α ωv (cid:88) k ∈R ( ¯ G ω ,v ) I ωk ( y ) + ¯ β ωv . (30)Obviously, the new Benders optimality cuts (29) and (30) are different from the old ones (27)and (28). able 10: Detailed statistics of instances MSG. Each row contains 5 instances.We report the number of solved instances (
S), the average CPU time in seconds(T), the average presolving time in seconds (PT), the average number of branchnodes (
N), and the average number of added Benders cuts (
C). TL denotesthat all five instances reach the time limit. (a)
ICM
Default SCNA INA K | Ω | p S T N C S T PT N C S T PT N C5 250 0.01 (b) LTM
Default SCNA INA K | Ω | S T N C S T PT N C S T PT N C5 250 able 11: Detailed statistics of instances GNU. Each row contains 5 instances.We report the number of solved instances (
S), the average CPU time in seconds(T), the average presolving time in seconds (PT), the average number of branchnodes (
N), and the average number of added Benders cuts (
C). TL denotesthat all five instances reach the time limit. (a)
ICM
Default SCNA INA K | Ω | p S T N C S T PT N C S T PT N C5 250 0.01 (b) LTM
Default SCNA INA K | Ω | S T N C S T PT N C S T PT N C5 250 able 12: Detailed statistics of instances HEP. Each row contains 5 instances.We report the number of solved instances (
S), the average CPU time in seconds(T), the average presolving time in seconds (PT), the average number of branchnodes (
N), and the average number of added Benders cuts (
C). TL denotesthat all five instances reach the time limit. (a)
ICM
Default SCNA INA K | Ω | p S T N C S T PT N C S T PT N C5 250 0.01 (b) LTM
Default SCNA INA K | Ω | S T N C S T PT N C S T PT N C5 250 able 13: Detailed statistics of instances ENRON. Each row contains 5 instances.We report the number of solved instances (
S), the average CPU time in seconds(T), the average presolving time in seconds (PT), the average number of branchnodes (
N), and the average number of added Benders cuts (
C). TL denotesthat all five instances reach the time limit. (a)
ICM
Default SCNA INA K | Ω | p S T N C S T PT N C S T PT N C5 250 0.01 (b) LTM
Default SCNA INA K | Ω | S T N C S T PT N C S T PT N C5 250 able 14: Detailed statistics of instances FACEBOOK. Each row contains 5instances. We report the number of solved instances (
S), the average CPU timein seconds (T), the average presolving time in seconds (PT), the average numberof branch nodes (
N), and the average number of added Benders cuts (
C). TLdenotes that all five instances reach the time limit. (a)
ICM
Default SCNA INA K | Ω | p S T N C S T PT N C S T PT N C5 250 0.01 (b) LTM
Default SCNA INA K | Ω | S T N C S T PT N C S T PT N C5 250 able 15: Detailed statistics of instances DEEZER. Each row contains 5 instances.We report the number of solved instances (
S), the average CPU time in seconds(T), the average presolving time in seconds (PT), the average number of branchnodes (
N), and the average number of added Benders cuts (
C). TL denotesthat all five instances reach the time limit. (a)
ICM
Default SCNA INA K | Ω | p S T N C S T PT N C S T PT N C5 250 0.01 (b) LTM
Default SCNA INA K | Ω | S T N C S T PT N C S T PT N C5 250 able 16: Detailed statistics of instances TWITTER. Each row contains 5instances. We report the number of solved instances (
S), the average CPU timein seconds (T), the average presolving time in seconds (PT), the average numberof branch nodes (
N), and the average number of added Benders cuts (
C). TLdenotes that all five instances reach the time limit. (a)
ICM
Default SCNA INA K | Ω | p S T N C S T PT N C S T PT N C5 250 0.01 (b) LTM
Default SCNA INA K | Ω | S T N C S T PT N C S T PT N C5 250 able 17: Detailed statistics of instances EPINIONS. Each row contains 5instances. We report the number of solved instances (
S), the average CPU timein seconds (T), the average presolving time in seconds (PT), the average numberof branch nodes (
N), and the average number of added Benders cuts (
C). TLdenotes that all five instances reach the time limit. (a)
ICM
Default SCNA INA K | Ω | p S T N C S T PT N C S T PT N C5 250 0.01 (b) LTM
Default SCNA INA K | Ω | S T N C S T PT N C S T PT N C5 250145.7 0 9320