Scalable and Memory-Efficient Algorithms for Controlling Networked Epidemic Processes Using Multiplicative Weights Update Method

1 Introduction

Due to the complex factors underlying disease spread, mathematical models based on the Susceptible-Infected-Recovered (or SIR) process are used extensively to understand the tradeoff between the cost of interventions (e.g., vaccination and social distancing), and the benefit (e.g., the number of infections averted), e.g., [Medlock and Galvani, 2009; Marathe and Vullikanti, 2013; Halloran et al., 2008; Eubank et al., 2004; Germann et al., 2006]. As an example, the CDC Scenario Hub [Truelove et al., 2021] involves the use of a variety of models-deterministic models based on differential equations, e.g., [Medlock and Galvani, 2009; Anderson and May, 1991], and stochastic models based on networks, referred to as networked SIR models, e.g., [Marathe and Vullikanti, 2013; Halloran et al., 2008; Eubank et al., 2004;

Germann et al., 2006; Chen et al., 2021]-in order to evaluate the benefits of different interventions, and find most effective ones.

We refer to the problem of designing an intervention to minimize the number of infections as EPICONTROL (note that this is a very general and complex problem, and we only study a very specific form of the problem, which is defined formally in Section 2). Such epidemic analyses implicitly try to find near-optimal solutions (or solutions better than current policies) to EPICONTROL to understand (a) how much benefit they can provide, or (b) structural properties of near-optimal solutions, which can give insights for more implementable policies, e.g., [Medlock and Galvani, 2009; Chen et al., 2021; Sambaturu et al., 2020]. We note that though the problem, in reality, is much more complex than the simplified version we study, even that is open, in general networks and disease model, making this a natural starting point.

Near-optimal solutions to EPICONTROL in models based on differential equations can be computed by local search methods, as in [Medlock and Galvani, 2009]. However, solving the EPICONTROL problem in network or agent based models is very hard [Eubank et al., 2006; Hayrapetyan et al., 2005]. A number of heuristics have been considered, e.g., prioritizing based on degree, centrality, or spectral properties e.g., [Cohen et al., 2003; Miller and Hyman, 2007; Saha et al., 2015; Chen et al., 2021], but it is possible to show that they have Ω(n) worst case bounds, in general, where n is the number of nodes in the network (we summarize other related work in Section 6). An approach that has been effective for some versions of the problem is using linear programming (LP) and randomized rounding [Hayrapetyan et al., 2005] (for the special case where the transmission probability is 1). This approach, combined with the sample average approximation (SAA) technique [Kleywegt et al., 2002] (which is a powerful tool from stochastic optimization, and reduces the problem to solving a deterministic problem on a set of samples) has been used for handling the more realistic regime of EPICONTROL with transmission probability less than 1 [Sambaturu et al., 2020; Babay et al., 2022]. [Babay et al., 2022] obtain the first rigorous approximation bounds for some settings; [Sambaturu et al., 2020] show that this approach works well in practice.

However, such LP based approaches (which use standard solvers, such as Gurobi) only scale to small county-sized networks (with up to 105 nodes). To put this in perspective, the CDC Scenario Hub models 1 run on state level network based models with millions of nodes. Further, such studies involve large experimental designs (due to the number of parameters), which requires solving intervention design problems on thousands of network models. Therefore, algorithms which can scale to networks with millions of nodes are necessary for supporting such policy analyses.

Key Contributions. We design algorithms MWUROUND and MWUROUND-SCALABLE for finding near-optimal vaccination strategies in networked SIR models in large networks with hundreds of millions of nodes (Section 3); see summary in Table 1. Our algorithms adapt the multiplicative weights update (MWU) [Arora et al., 2012] to run the sample average approximation and LP rounding steps of [Sambaturu et al., 2020]. More specifically,

• We show that our problem can be reduced to a covering type problem, for which prior MWU based algorithms exist [Arora et al., 2012; Fleischer, 2000]; we show that by exploiting the problem structure, we can get an improvement of more than a factor of M (the number of samples) over the standard use of MWU (Section 3). This also improves both the running time and space of SAAROUND by more than a factor of mM-this is many orders of magnitude in massive networks.
• We design a more scalable algorithm MWUROUNDSCALABLE, by running the MWU computations on random samples computed on-the-fly, leading to almost a factor M improvement in running time, and a factor M improvement in space. This allows MWUROUNDSCALABLE to run on national scale networks.
• We show that our approach can be extended to the FAIREPICONTROL problem, which incorporates fairness constraints with respect to the budget. Our results show that such fairness constraints have a significant impact on epidemic outcomes.

These performance gains (runtime and memory) are achieved without compromising on the approximation quality. In Section 5, we present a detailed experimental evaluation of the proposed algorithms on a number of real-world and synthetic networks; the largest is a country-scale contact network containing over 334M nodes. Though the version of EPICONTROL we consider is admittedly the simplest possible (which itself has not been fully resolved), e.g., compared with other models, our results show that the sample average approximation, linear programming and rounding, and MWU techniques can help in scaling interventions in networked SIR models.

1 https://viz.covid19scenariomodelinghub.org

The supplementary information, including proofs and further details, can be found in the full version 2 .

2 Preliminaries

Let G = (V,E) be a contact network, where V is the set of people or nodes, and e = (u, v) ∈ E when nodes u, v ∈ V come into direct contact. Let |V | = n and |E| = m. Let us assume the following simple SIR model of disease spread.

SIR Epidemic Model on Networks. An agent based SIR model of disease spread can be viewed as a diffusion process on a network [Marathe and Vullikanti, 2013; Halloran et al., 2008]. Each node in the network is in one of the following three states: Susceptible (S), Infectious (I), or Recovered/Removed (R). In its simplest form, we have a discrete time model, and in each time step, an infected node spreads infection to each of its susceptible neighbors with a probability p - called transmission probability. An infected node recovers after its infectious duration (and doesn't get infected after that). We assume that the disease starts at a fixed set of externally infected nodes S ⊆ V and |S| = k (our results extend to more general starting conditions, e.g., a random source). We assume that interventions have 100% efficacy. Let #infections(X) denote the number of infections that result if set X of nodes is vaccinated; this is a random variable, and our goal is to minimize E[#infections(X)].

EPICONTROL problem. Given a contact network G = (V,E), a simple SIR model of disease spread, a fixed set of sources of infection S where |S| = k, and a budget B on the number of interventions. Find a set of nodes X ⊆ V such that |X| ≤ B and the expected number of infections E[#infections(X)] is minimized.

Our results extend to more general settings with nonuniform transmission probabilities, and other starting conditions, such as random sources. However, the EPICONTROL stated above is admittedly much simpler than the ones that need to be considered in public health analyses, e.g., requiring adaptive interventions, which have limited efficacy, and with a lot of uncertainty in the model.

2 Link to full version: https://tinyurl.com/mwuroundfullversion

Next, we consider fairness with respect to the budget. We have a partition of V into groups V1, . . . , Vc, and a budget Bi for each group (e.g.,proportional to the size of each Vi).

FAIREPICONTROL problem. Given a contact network G = (V,E), a simple SIR model of disease spread, a fixed set of sources of infection S where |S| = k, a partition of V = ∪ci=1Vi, and a budget Bi for each group Vi. Find a set of nodes Xi ⊆ Vi such that |Xi| ≤ Bi, for each i, and the expected number of infections E[#infections(X1, . . . ,Xc)] is minimized.

Note that there are many other more complex fairness models [Barocas et al., 2019], but the EPICONTROL is challenging even in the simplest setting, e.g., [Dinitz et al., 2022].

2.1 SAA Based Algorithm

We briefly summarize the stochastic optimization based algorithm of [Sambaturu et al., 2020], referred to as SAAROUND, which is our starting point. It involves the following steps (we only give a short description here, and refer to [Sambaturu et al., 2020] for details).

• Construct a sampled graph H′j = (VH′j ,EH′j ), for j = 1, . . . , M, by picking each edge e ∈ E to be in EH′j with a probability p.
• Solve the following linear program (LPsaa): min 1M ∑j ∑v yvj, subject to: (1) ∀j, ∀u ∈ V : yuj ≤ 1-xu, (2) ∀j, ∀u ∈ V, (w, u) ∈ Ej : yuj ≥ ywj -xu, (3) ∀s ∈ S : ysj = 1, (4) ∑u∈V xu ≤ B (5) All variables ∈ [0, 1]
• Let x, y be the optimal fractional solution to LPsaa. For each v, set Xv = 1 with probability min{1, 2xv log(4nMN)}, where N is the maximum number of paths from S to any node v in H′j .
• X = {v : Xv = 1} is the set of nodes vaccinated.

[Sambaturu et al., 2020] show that using M = Ω(n2 log n) samples, the above approach gives a bicrteria approximation, vaccinating O(log(4nMN)B), while ensuring that the expected number of infections is at most six times the optimal. Solving the LP in SAAROUND is the main bottleneck, as it has n+nM variables and ∑j |Ej| constraints-they report that this only scales to contact networks with up to 105 nodes.

3 Algorithm

We improve the approach of [Sambaturu et al., 2020] by bypassing the need to use a solver to directly solve the LP. Instead, we adapt the multiplicative weights update technique [Arora et al., 2012] to find a near-optimal solution to the LP. It will be easier to present the LP in a slightly different form. Let yvj be an indicator whether the node v gets infected in the sampled outcome H′j . Let P ′v be the set of paths from S

to node v in any outcome H′j for j ∈ [M].

By adding up the constraints for each edge of a path of LPsaa, it can be verified that we get the constraint (2) of LPpath, which is summarized below.

Observation 1. The above LP is equivalent to LPsaa.

The main ideas and steps are described below.

1. Lagrangian multiplier for budget. The dual of LPpath is complicated due to a negative coefficient associated with the budget constraint (3). We simplify it by changing the objective to 1M yvj + λ ∑u∈V xu, with the multi-
2. ∑j ∑v∈VH′j

plier λ for the cost of the solution. The budget constraint is dropped; we refer to this LP as LPLM. This simplifies the resulting LP, since it only has covering constraints. As λ increases, ∑u xu will decrease in the optimal solution. Since values for λ are not known a priori, a binary search can be employed to find a suitable value λ′ such that ∑u xu ≤ B, which is done in Algorithm LSEARCH-SAA. The x′ , y′ values returned by Algorithm 1 for λ′ provides an approximate solution to LPsaa.

1. Constructing augmented sampled graphs. For simplifying the presentation, we construct M sampled graphs Hj = (VHj ,EHj ) in the following manner: Hj is initially the same as H′j , constructed as in the first step of SAAROUND. Let A(j) = {a(u, j) : u ∈ VH′j - S}, where a(u, j) denotes a copy of node u in Hj, and is referred to as a stub. Let A = ⋃j∈[M] A(j) be set of all stubs. Each stub a(u, j) is attached to u by an edge (u, a(u, j)). Overloading the definitions, let Pv,j denote the set of paths from S to a stub node v = a(u, j) ∈ A(j) in Hj. Let Pj = ⋃v∈A(j) Pv,j and P = ⋃j∈[M] Pj .
2. Variables and costs. We associate a length to each node u ∈ (V - S) ⋃ A denoted by ℓ(u). For a node u ∈ V , ℓ(u) will correspond to the variable xu, while for a node v = a(u, j) ∈ A, ℓ(v) will correspond to the variable yuj. The length of any path in P ∈ P is given by the sum of lengths of nodes on this path. Let ℓ = ⟨ℓ(u) : u ∈ V ⋃ A⟩ denote the vector of length variables. Let c(u) = λ for u ∈ V \S denote its capacity, whereas c(u) = 0 for u ∈ S. Let c(v) = 1M for v = a(u, j) ∈ A denote the capacity of a stub v. We will keep track of flows on the network; let z(P) denote the flow on the path P ∈ P. Let z = ⟨z(P) : P ∈ P⟩ denote the vector of flow variables.

Algorithm 1 MWUSAA (λ)

Input: parameter λ (we assume the network G = (V,E), S, subgraphs H1, . . . ,HM, ϵ are fixed, δ = (1 + ϵ)((1 + ϵ)L)- 1ϵ where L is the max. number of nodes on any path in G Output: ℓ

1: Initialize ℓ(u) = δ for all u ∈ (V -S) ∪A, z(P) = 0 for all P ∈ P. 2: Set c(u) = λ for u ∈ V - S and c(v) = 1/M for v ∈ A 3: for r = 1 to ⌊log1+ϵ 1+ϵ δ ⌋ do 4: for j = 1 to M do 5: while there exists pathP ∈ Pj such that ℓ(P) < δ(1+ϵ)r do 6: Let c(P) = minu∈P c(u) 7: Let d ≥ 1 be the smallest integer such that ∑v∈P-S ℓ(v)(1 + ϵc(P) c(v) )d ≥ δ(1 + ϵ)r 8: z(P) ← z(P) + d · c(P) 9: For v ∈ P - S, ℓ(v) ← ℓ(v)(1 + ϵc(P) c(v) )d 10: end while 11: end for 12: end for 13: for each v ∈ V - S, ℓ(v) = ℓ(v) ℓmax where ℓmax = maxv∈V \S ℓ(v) 14: Return ℓ

1. Simplified LP. The above discussion reduces our LP to the following:

5. Incremental computation of ℓ(·). AlgorithmMWUSAA computes an approximate solution to LPℓ, using the multiplicative weight update technique [Arora et al., 2012]. It starts by initializing the length ℓ(v) = δ for each u ∈ (V \ S) ⋃ A, where δ has a very small value determined in the analysis. The ℓ(v) for v ∈ S is initialized to zero. Also, for each u ∈ V - S, we set a capacity c(v) = λ, whereas for v ∈ A the capacity c(v) = 1M . In each iteration r of the algorithm, and for each augmented sampled graph Hj, we update the lengths of nodes on the paths in Pj corresponding to an augmented sampled graph Hj until all the paths in Pj have length at least δ(1 + ϵ)r - this value is referred to as threshold(r) for the rth iteration. The algorithm terminates after rmax = ⌊log1+ϵ 1+ϵ δ ⌋ iterations. Since, the threshold(rmax) for the rmax iteration is in [log1+ϵ 1+ϵ δ - 1, log1+ϵ 1+ϵ δ ], we are guaranteed that, at termination, all paths in P are of length in range [1, 1 + ϵ], thereby satisfying the constraints of the linear program.

Subroutine. LSEARCH-SAA(M,B) starts with λ = 1 MB , uses binary search to find largest λ that has ∑u∈V \S ℓ(u) ≤ B; compute ℓ = MWUSAA(λ), and increase λ if the condition on ∑u ℓ(u) is still satisfied (see full version for details).

Algorithm 2 MWUROUND(G, S, M, B, p, ϵ)

• 1: ℓ = LSEARCH-SAA(M, B)
• 2: Using the randomized rounding in [Sambaturu et al., 2020], round the fractional solution ℓ to an integral solution X
• 3: X = {u : u ∈ V \ S and X(u) = 1} is the set of nodes picked for intervention
• 4: return X

Analysis. We show below that Algorithm MWUSAA gives an approximate solution to LPℓ(λ). Theorem 1 summarizes its running time. Note that this is a factor M better than directly using the MWU technique of [Fleischer, 2000].

Theorem 1. MWUSAA can be implemented using O(n log1+ϵ 1+ϵ δ ) shortest path computations, and has a total running time of ˜O(ϵ-2nmM).

3.1 Incorporating Fairness

To adapt our approach to solve FAIREPICONTROL, we have a separate lagrangian multiplier λi for group Vi, and the objective becomes 1M ∑j ∑v∈VH′j yvj + ∑ci=1 λi ∑u∈Vi xu. As

in Algorithm MWUSAA, we search over the space of the λi's, till the budgets are all within the required bounds.

4 AlgorithmMWUROUND-SCALABLE

The main bottleneck in MWUSAA is that in each phase, it has to iterate over all the M sampled graphs, thereby needing to stores all the M graphs in memory. This can be handled by storing the sampled graphs in files and loading only one graph into memory at any time, however, this affects the total runtime. Our next algorithm MWUSCALABLE, whose main ideas are presented below, addresses these issues.

1. Generate random stubs. The intuition behind this approach is that the actual samples do not matter, as long as we are able to generate the paths that would appear in these samples in each iteration of the algorithm. Let the probability that u is reachable from S in a sampled graph be denoted by sp(u), and is referred to as stub probability. This can be estimated from our sampling process as follows: sp(u) ≈ reachable(u,S,M) M , where reachable(u, S,M) denotes the number of samples in the M sampled graphs in which u is reachable from sources S. At the start of MWUSCALABLE, we generate the random set of stubs Asp(u) for each node u as follows: for each u ∈ V and j ∈ 1, · · · ,M, the stub a(u, j) is generated with probability sp(u). Let Asp = ⋃u Asp(u). The initial length ℓ(v) = δ for each v ∈ Asp.

1. Generate sampled graphs on the fly. In every iteration r, we generate only q ≪ M sampled graphs, H′j = (VH′j , EH′j ) for j ∈ [1, q]. The algorithm then works on one sampled graph at a time, therefore, at any time stores only one sampled graph in memory.
2. Phases and iterations of the algorithm. In each phase, q iterations are performed. In each iteration q of phase r, the algorithm generates a random sampled graph H′j . Then, independently attaches a random stub for each node reachable from sources in H′j to form Hj. Then, it iteratively updates lengths of paths in the sampled graph Hj until there is no path of length less than threshold(r). The algorithm terminates after rmax = ⌊log1+ϵ 1+ϵ δ ⌋ phases (same as MWUSAA). The pseudocode of this algorithm is provided in the full version.
3. Computation of y values. Some of the ℓ(u) for u ∈ V could have a value in [1, 1 + ϵ]. Therefore, to make the solution x feasible, we scale ℓ(u) = ℓ(u) ℓmax for u ∈ V where ℓmax = maxu∈V ℓ(u). Since the sampled graphs Hr generated in each iteration are a combination of paths from many sampled graphs, the ℓ(v) variables for v ∈ Asp will not be meaningful. Therefore, we use the ℓ(u) for u ∈ V obtained after scaling, and re-compute ℓ(v) (corresponding to yvj variables) for v = a(u, j) as follows: for each sampled graph H′j , find a shortest path tree of H′j with S as sources using ℓ(u) for u ∈ V \ S as weights. For each u ∈ VH′j , let Puj be a shortest path to node u that has length ℓ(Puj). Then, for the stub v = a(u, j), we set ℓ(v) = min{0, 1 - ℓ(Puj )}.

Lemma 1. The solution ℓ computed by MWUSCALABLE is a feasible solution to LPℓ(λ).

The search over λ values is an embarrassingly parallel task, as these searches are independent of each other. Therefore, we implement a parallel version of this algorithm.

5 Experiments

We address the following questions in our experiments:

1. Performance. What are the empirical guarantees of our methods? How does the performance of our approach compare to the baselines for this problem?
2. Impact of parameters. How do the runtime and solution quality of our methods vary with changes in transmission probability p and error parameter ϵ?
3. Scaling. How does the runtime of our approach grow with that of size of the network?
4. Parallelism. What is the throughput of our parallel implementation?
5. Cost of fairness. What is the cost of incorporating fairness constraints?

ϵ

Figure 2: Comparison of approximation ratios of fractional solutions obtained by MWUROUND and MWUROUND-SCALABLE. B = 50.

5.1 Datasets andMethods

In our experiments, we considered networks of different classes and varying sizes (Table 3) for evaluating the performance and scalability of our approach-these range from random networks (PA1) based on the preferential attachment model [Barab´asi and Albert, 1999], collaboration networks, such as CA-GrQc and CA-HepTh [Leskovec et al., 2007]), synthetic contact networks for Montgomery county VA and Portland city from [Sambaturu et al., 2020], and for Virginia from [Chen et al., 2021], and networks Regional and US-size - which are generated by connecting many copies (5 and 44 respectively) of the Virginia network using random edges, in order to study scaling.

Methods and Baselines. We consider the following methods: (1) SAAROUND [Sambaturu et al., 2020], (2) DEGREE: select top B nodes in V \S for intervention, (3) NO-ACTION: Baseline with no interventions, (4) MWUROUND, and (5) MWUROUND-SCALABLE and its parallel implementation.

5.2 Performance

The full version provides additional details and plots.

Figure 2 shows that the approximation ratio of fractional solutions obtained by MWUROUND is within 1.2, even for ϵ = 0.15. In comparison, the approximation ratio obtained by MWUROUND-SCALABLE (for q = 1) is at most 1.3 for ϵ = 0.04, and goes up to 1.7 for ϵ = 0.15, which is within a factor of (1 + 5ϵ). MWUROUND-SCALABLE has significantly better performance for small values of ϵ and on higher values of p (full version).

5.3 Runtime Performance

Figure 3 presents the runtime performance (for a fixed λ) of MWUROUND-SCALABLE on various networks. The number of λ values considered by MWUROUND-SCALABLE algorithm determines the runtime of the algorithm. The depth of λ search decreases with the increase in the budget B. MWUROUND-SCALABLE ran within a few minutes (< 15 minutes averaged over a few runs), for a fixed λ value, on the Portland network for problem, whereas it ran in just about 2 days on the US-size network which has over 334 million nodes and 32 billion edges - for instances with a high attack rate (> 40% population infected).

5.4 Parallel Implementation

Figure 4 shows how the throughput of MWUROUNDSCALABLE changes when increasing the number of threads and varying the budget on Virginia. The peak throughput observed, is between 16 and 32 threads, for our computing platform. At its maximum throughput, MWUROUNDSCALABLE scales graciously with the size of the network.

5.5 Cost of Incorporating Fairness

The nodes in Montgomery network are divided into two age groups: (i) nodes with age at most 19, and (ii) nodes with age over 19. We consider two settings: (i) without fairness: budget B is allotted to all nodes, and (ii) with fairness: budget B is allotted to each group as a proportion of its size. Figure 5 shows the ratio of the expected number of infections in a fair allocation to that without fairness-observe that this ratio can be as high as 1.4, which is the price of fairness.

5.6 Discussion

MWUROUND-SCALABLE outperforms the degree baseline, and is faster than the SAAROUND which uses the LP solver. These runtimes can be improved using our parallel implementation. Both, MWUROUND-SCALABLE and its parallel version, scale well for networks larger than Portland, and run even on the US-size network which has over 32 billion edges.

6 RelatedWork

There is a huge amount of literature on interventions for epidemic models. Due to the limited space, we only discuss network based models [Marathe and Vullikanti, 2013; Halloran et al., 2008; Lofgren et al., 2014; Eubank et al., 2004; Germann et al., 2006], which is our focus here; see full version for further details. Such models have been found to be more powerful and useful for epidemic spread on large heterogeneous populations, where the complete mixing assumptions of differential equation models are not valid. However, these are harder to set up, simulate and optimize over. Therefore, a number of heuristics have been proposed, which prioritize vaccination based on degree, centrality, or spectral properties e.g., [Cohen et al., 2003; Miller and Hyman, 2007; Saha et al., 2015; Minutoli et al., 2020]. These do not directly give rigorous worst case guarantees; however, some of them can be computed very efficiently. The influence based approach of [Minutoli et al., 2020] has been parallelized using clever hill climbing techniques.

7 Conclusions

Our algorithms are the first to give near-optimal vaccination strategies and minimize the expected number of infections in SIR models on US-scale networks. The reduced time and memory requirements result from the multiplicative weights update method adapted for percolation processesin contrast, prior methods, which used the state-of-the-art LP solvers, were not able to scale beyond county-scale contact networks. Further improving the performance for US-scale networks is a natural question. We have focused on nonadaptive interventions; extending this to more general interventions is a challenging open problem. Finally, uncertainty is an important component of all such models. The SAA technique is easily amenable to handling some models of uncertainty. Extending our work to incorporate uncertainty is another important problem.

Acknowledgments

We would like to thank all the funding agencies for supporting this work. This work has been partially supported by the following grants: DTRA CNIMS (Contract HDTRA111-D-0016-0001), NSF DIBBS Grant ACI-1443054, NSF EAGER Grant CMMI-1745207, NSF BIG DATA Grant IIS-1633028, NSF Expeditions in Computing Grant CCF1918656, NIH R01GM109718, and the Exascale Computing Project's (ECP)) (17-SC-20-SC) ExaGraph codesign center. ECP is a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration. The research was in part supported by the NSF awards OAC 1910213 and CCF 1919122, and USDA NIFA award 2021-67021-35344.