Fast Combinatorial Algorithm for Optimizing the Spread of Cascades

1 Introduction

Spatial conservation planning problems have received significant attention from the AI community, resulting in a range of strategies for conserving land parcels in order to support the recovery of an endangered species or preserve biodiversity. Sheldon et al. [2010] studied a restricted version of this problem-Red-cockaded Woodpecker (RCW)-in which the planner selects a set of land parcels, subject to a budget constraint, to maximize the spread of a population over a geographical network of land patches within a specific time horizon. The underlying spreading process is modeled as a network diffusion process or cascade [Kempe et al., 2003] using a metapopulation model developed by ecologists. Due to the stochasticity of process, the RCW problem can be written as a stochastic optimization problem. Similar propagation models have been studied in social networks, particularly the influence maximization problem where a planner strategically selects a certain number of sources to trigger the spread of influence over a social network [Domingos and Richardson, 2001; Kempe et al., 2003]. Unlike that problem, a major challenge presented by RCW is that its objective function is not submodular. Consequently, greedy algorithms may produce solutions that are arbitrarily worse than an optimal solution.

Sheldon et al. [2010] developed an algorithm based on sample average approximation (SAA) and mixed integer program (MIP) to solve the problem approximately. In SAA, a certain number of random samples are drawn and their average is used to approximate the original objective function that is an expectation over all possible cascade scenarios. This way, the RCW problem is converted into a discrete optimization problem and solved optimally by a standard MIP solver. However, since the discrete optimization problem is NP-hard, the MIP solver can only handle limited network sizes and a small number of samples, undermining the approximation quality of SAA. For example, a standard MIP solver fails to work with more than 20 samples, either taking too long to finish or using too much memory, and an improved MIP based algorithm [Kumar et al., 2012] can scale up to 40 samples but needs hours to finish. To address this challenge of scalability, we propose what we believe is the first combinatorial algorithm to solve the converted discrete optimization problem approximately and efficiently.

In fact, our algorithm is applicable to a more general problem called Budget Set Weighted Directed Steiner Graph (BSW-DSG) and, with minor modifications, to the Quota Set Weighted Directed Steiner Graph (QSW-DSG) and PrizeCollecting Set Weighted Directed Steiner Graph (PCSW-DSG) problems. We introduce each of these problem as a generalization of a classical network design problem to the setting where edges are purchased in sets instead of individually. Xue et al. [2015] first generalized the Steiner tree problem in this way; our work extends this generalization to the widely-used quota, budget, and prize-collecting Steiner tree variants [Johnson et al., 2000]. All three problems share the same basic setup but have different constraints and objectives. In the basic setup, an initial graph and a set of candidate edge sets, each with a certain cost, are given. The planner can purchase a subset of edge sets and add the included edges into the initial graph. In the budget problem (BSW-DSG), the total cost of selected edge sets is constrained by a given budget and the objective is to maximize the number of terminal vertices connected to a unique root vertex in the augmented graph. In the quota problem (QSW-DSG), the number of terminal vertices connected to the root is required to be greater than a given value and the objective is to minimize the total cost of achieving this requirement. In the prize-collecting problem (PCSW-DSG), both the number of terminal vertices connected to the root and cost are part of the objective and the tradeoff between them is controlled by a parameter. In this paper, we propose a fast primal-dual algorithm to solve the PCSWDSG problem approximately. To address the BSW-DSG problem, we use Lagrangian relaxation [Jain and Vazirani, 2001; Ahuja et al., 1988] to fold the budget constraint into the objective and convert the BSW-DSG problem into a PCSW-DSG problem that is parameterized by a Lagrangian multiplier β. We also derive a bisection procedure to find the value of β that gives a near optimal solution. At each iteration of the procedure, a PCSW-DSG problem parameterized by some β is solved near optimally. Similarly, the same algorithm can solve the QSW-DSG problem with a minor modification to the bisection procedure. For the sake of clarity, we focus the paper on the BSW-DSG problem.

When applied to the RCW dataset, our algorithm produces near optimal solutions and runs orders of magnitude faster than the MIP approach, for example, taking only 20 minutes for 300 samples. A greedy baseline also produces near-optimal results, however our algorithm still produces better quality in most cases and, most importantly, is 10-30 times faster. In experiments on synthetic problems with objective functions that are strongly non-submodular, our algorithm performs much better than the greedy baseline.

The rest of the paper is organized as follows. Section 2 describes the conservation planning problem and the basic solution based on SAA. Section 3 formalizes the general BSWDSG problem. Section 4 presents our combinatorial algorithm. Section 5 describes how to apply our algorithm to a range of other practical applications. Experimental results are reported in Section 6.

2 Problem Statement

Consider a conservation area consisting of habitat patches that are the atomic units in the population dynamics model and can be either occupied or unoccupied by the species. These patches are grouped into non-overlapping parcels P1, P2, ..., PL, which are the smallest units available for purchase. A parcel Pl can be purchased with cost cl used to restore and conserve its habitat patches so they are suitable for the species to occupy. A patch can be occupied only if the parcel that it belongs to is purchased. A conservation strategy is an L-dimensional vector y where yl = 1 if that the lth parcel is purchased and yl = 0 otherwise.

Sheldon et al. [2010] use a metapopulation model from ecology [Hanski and Ovaskainen, 2000] to describe the stochastic occupancy dynamics of the species. In this model, a patch is either occupied or unoccupied at each time step h ∈ {0, 1, ...,H - 1}. At time 0, the species occupy a set of patches called sources. At any time h, an occupied patch v triggers two stochastic events: (1) the population at u colonizes an unoccupied patch v with probability puv (typically decays with distance between the patches) that makes v occupied at time h+1; and (2) the population at v becomes extinct with probability 1 - pvv, making v unoccupied at time h + 1.

The goal of conservation planning is to select a set of parcels to purchase, without exceeding a budget B, such that the expected number of patches being occupied at time H-1 is maximized. More concretely, given a strategy y, define 0-1 random variable X(y)v,h for each patch v and each time h to indicate whether the patch v is occupied at time h (X(y)v,h = 1) or not (X(y)v,h = 0). The best strategy is obtained by solving the following optimization problem:

Sample Average Approximation The stochastic optimization problem (1) is very hard to solve directly, but the SAA scheme can be used to solve it approximately. The idea is to sample N independent cascade scenarios from the stochastic spreading process, each of which represents a possible outcome of the stochastic process. For each scenario, given a strategy vector y, the occupancy status of any patch at any time step can be evaluated directly so that the number of occupied patches at the last time step can be determined. Hence, the average number of occupied patches over sampled scenarios is an approximation of the objective in (1).

More concretely, one scenario is a layered graph, as shown in Fig. 1, that contains a vertex vh for patch v and time h. To construct the edges of the graph, for a pair of patches (u, v) (u ̸= v) and a h ∈ {0, . . . , H - 2}, a biased coin with probability puv being heads is flipped. A directed edge (uh, vh+1) is added to the graph iff the outcome is heads, meaning that v will be occupied at h + 1 if u is occupied at h. Similarly, for a patch v and h ∈ {0, . . . , H - 2}, a coin with probability 1 - pvv being heads is flipped. A directed edge (vh, vh+1) is added to the graph iff the outcome is heads, meaning that the population occupying v at h will survive to time h+1. With a conservation strategy, we can determine whether v is occupied at h by checking whether there exists a valid path from some u0 to vh such that (1) u is occupied at time 0 and (2) all patches on the path are purchased. By checking the layered graphs of N i.i.d. sampled cascade scenarios, we can obtain the average number of occupied patches at time H - 1. Specifically, let X(y)nv,h be a 0-1 variable indicating whether the patch v is occupied in the nth scenario at time h (X(y)nv,h = 1) or not (X(y)nv,h = 0). Then, we have

When N goes to infinity, the right term converges to the left term and the vector y that maximizes the right term converges to the vector that maximizes the left term [Kleywegt et al., 2002]. Thus, a good conservation strategy can be obtained by solving the following optimization problem:

Instead of solving (3) as a MIP, we propose a significantly faster near-optimal algorithm.

3 Budget Set Weighted Directed Steiner Graph Problem

In this section, we propose a novel network design problem called Budget Set Weighted Directed Steiner Graph (BSWDSG), of which (3) is a specific instance.

The input of a BSW-DSG problem is a directed graph G = {V,E} with a unique root vertex r ∈ V , a set of terminal vertices T ⊆ V and a collection of M edge sets E = {E1, E2, ..., EM} where Es ⊆ E and each Es is associated with a nonnegative cost cs. Let A denote a subset of E. Each A corresponds to a subgraph GA = {V, EA} where EA = ⋃Es∈A Es. A vertex v is connected to r if there is a path from r to v. Given a budget B, the goal is to find a set A with cost no greater than B such that the number of terminal vertices connected to r in GA is maximized.

To write down its mathematical formulation, we start with a preprocessing step. Let ET be an additional set of edges from the root directly to each terminal k ∈ T . For any selected set A, we can augment the subgraph GA to be a Steiner graph (one where all terminals are connected to r) by adding edges (r, k) ∈ ET for each terminal k that is not already connected to r in GA. Since T is fixed, the original goal of maximizing the number of connected terminals is equivalent to minimizing the number of edges needed to augment GA to obtain a Steiner graph. Mathematically, BSW-DSG is formulated as follows:

where AT is the set used for augmentation.

Problem (4) can be formulated as a MIP using the networkflow encoding shown in Fig. 2. The flow variable xku,v encodes the flow destined for terminal k ∈ T on the edge (u, v) ∈ E, where E denotes E ⋃ ET . It is easy to show that, in the optimal solution of the MIP, xku,v is either 0 or 1. The 01 decision variable ys indicates whether the edge set Es is purchased (ys = 1) or not (ys = 0). The 0-1 variable zrk indicates whether the edge (r, k) ∈ ET is used (zrk = 1) or not (zrk = 0). Line (2) forces each terminal k to be connected to r by requiring one unit of flow on a path from r to terminal k. The objective is to minimize the number of used edges in ET . Line (3) says that the flow can pass an edge (u, v) only

Figure 2: MIP formulation of BSW-DSG problem

when (u, v) is in a purchased edge set Es, (4) says that the flow can pass an edge (r, k) ∈ ET only when rrk is 1, and (5) is the budget constraint.

BSW-DSG Formulation of the RCW problem The optimization problem (3) can be formulated as a BSW-DSG problem in the following way. A new input graph of BSW-DSG is created by combining the layered graphs of N sampled scenarios in which a vertex vn h exists for each patch v, each scenario n and each time step h. Also, a new root vertex r is added to the graph along with directed edges from r to each vn 0 with v being occupied at time 0. Then, vn h being connected to r implies the patch v is occupied at time h in the nth scenario. The terminal set T is defined to contain all vertices with subscript H - 1, that is, T = {vnH-1 : ∀v, n ∈ {1, ...,N}} since the goal is to optimize occupancy of patches at time H - 1. An edge set Es is created for a parcel Pl in a way that if a patch v is in Pl, all edges ending with vn h is in Es. Formally, Es =

The cost of Es is the cost of Pl. By now, we formulated the optimization problem (3) of the RCW problem as an instance of BSW-DSG.

4 Our Algorithm for BSW-DSG Problem

Lagrangian relaxation has been used to solve many constrained discrete optimization problems [Jain and Vazirani, 2001; Kumar et al., 2012; Ahuja et al., 1988]. The basic idea is to relax a complex constraint and bring it to the objective together with a Lagrangian multiplier β. The new optimization problem is called relaxation problem parameterized by β. In this paper, we take advantage of the problem structure to develop a bisection procedure to search for a β such that the near optimal solution can be computed by solving the relaxation problem with that β. For the BSW-DSG problem, the relaxation problem is a PCSW-DSG problem with tradeoff parameter β. To solve the PCSW-DSG problem, we propose a primal-dual algorithm motivated by [Wong, 1984].

In summary, our algorithm consists of three major steps. First, we relax the budget constraint (5) in Fig. 2 and move

PCSW-DSG Problem / Primal Problem:

Figure 3: Primal and dual formulations of relaxation problem

it to the objective with a parameter β to create the PCSWDSG problem shown in Fig. 3. Second, we use the bisection procedure to find the β, using the primal-dual algorithm to solve the PCSW-DSG problem. Finally, we extract and refine the solution. The primal-dual algorithm and the bisection procedure will be explained next.

4.1 Bisection Procedure

Let's first analyze the function L(x, y, z, β) in Fig. 3. Let (xβ, yβ, zβ) denote the optimal solution of the PCSW-DSG problem for some β. Let, Z(β) = ∑(r,k′ )∈ET zrk′ β, C(β) = ∑Ms=1 csysβ, and L(β) = Z(β) + βC(β). We have the following properties.

Proposition 1. As β increases, Z(β) is nondecreasing and C(β) is nonincreasing.

The proof is omitted due to space limitation. Intuitively, larger β puts larger penalty on the cost and thus results in a less costly y but larger Z(β).

The bisection procedure starts with a suitably large interval bounding β and then narrows it iteratively. At each iteration, a new β is picked as the middle point of the current range and then the primal-dual algorithm is used to calculate C(β) and produce a near optimal strategy for β. By Proposition 1, if C(β) < B, all β′ greater than β will give Z(β′ ) ≥ Z(β) and therefore can be abandoned. If C(β) > B, all β′ smaller than β will give C(β′ ) ≥ C(β) and therefore can be abandoned. The procedure terminates when the range is less than some threshold and then a β in the range along with the computed strategy is returned.

Proposition 1 is true when the PCSW-DSG problem is solved optimally for each β, but the bisection procedure is still valid with our approximate algorithm. Let (ˆxβ, ˆ yβ, ˆ zβ ) denote the near optimal solution computed by the primal-dual algorithm for some β. And let ˆ Z(β) = ∑(r,k′ )∈ET ˆ zrk′ β and ˆ C(β) = ∑Ms=1 cs ˆ ysβ. Then, we have

Proposition 2. Suppose terminal vertices in line 6 of Algorithm 1 are chosen in a predefined order, as β increases, ˆ Z(β)

Algorithm 1 Primal-Dual Algorithm for PCSW-DSG

1: function PRIMALDUAL(G, E, β) 2: λkuv ← 0 ∀k ∈ T ∀(u, v) ∈ E 3: G′ = {V,E′ }, E′ ←φ 4: y ← 0, z ← 0 5: while G′ is not a steiner graph do 6: Pick a k ∈ T not connected to r 7: Find δ(k) : a set of cut edges between r and k 8: S = {s|Es ∩ δ(k) ̸= 0, Es /∈ A} 9: s∗ = arg mins∈S ∆(s, k) where 10: ∆(s, k) = βcs - ∑k′∈T ,(u,v)∈Es λk′ uv/|Es ∩ δ(k)| 11: if ∆(s∗ , k) ≥ (1 - λkrk) then 12: ∆(s∗ , k) ← 1 - λkrk 13: E′ ← E′ ∪ {(r, k)} and zr,k ← 1 14: else 15: E′ ←E′ ∪Es∗ and ys∗ ← 1 16: end if 17: λkuv ← λkuv +∆(s∗ , k) ∀(u, v) ∈ δ(k) 18: end while 19: Set ys ← 0 if removing Es, G′ is still Steiner graph. 20: Return y 21: end function

is nondecreasing and ˆ C(β) is nonincreasing.

The proof is omitted. Intuitively, as β becomes larger, the equality in (5) in the dual problem becomes harder to satisfy. Since the terminal vertices are picked in the same order over different β, the set of edge sets purchased for a larger β is the subset of edge sets for a smaller β which gives us Proposition 2. Finally, if the primal-dual algorithm performs near optimally for the PCSW-DSG problem, the final solution computed by the bisection procedure will be near optimal for the BSW-DSG problem. Although this procedure is not accompanied by a formal approximation guarantee, we observe that it produces near-optimal results on the problem instances used in our experiments.

4.2 Primal-Dual Algorithm for PCSW-DSG

Algorithm 1 summarizes the PCSW-DSG solver. This algorithm builds technically on the primal-dual algorithm of Xue et al. [2015] for the SW-DSG problem (where all terminals must be connected). The dual formulation, shown in Fig. 3, is obtained based on the LP relaxation of the primal problem. The basic idea of the primal-dual algorithm is to repeatedly adjust dual variables to increase the dual objective and at the same time construct a feasible primal solution while maintaining the primal complementary slackness condition [Vazirani, 2001]. This condition is defined as follows for our problem.

Definition 1. Primal Complementary Slackness Condition:

Without loss of generality, we assume that E doesn't contain edge (r, k) for any k ∈ T . Initially, the primal variables y, z are set to be 0 (line 4), meaning that no edge sets are purchased and no edges in ET are used (line 3), which represents an infeasible solution for the primal problem. Also, the dual variables are set to be 0 (line 2), which is a feasible solution for the dual problem. Note that we omit the dual variables {µkv } in the presentation of the algorithm because their values can be derived from the values of {λkuv } and our goal is to construct the feasible primal solution rather than explicitly calculate the actual value of the dual objective. The relationship of µ and λ is described by the following Proposition.

Proposition 3. Consider the edge-weighted graph Gk = {V, E} where each edge (u, v) is assigned length equal to λkuv . For any pair of vertices u and v where there is a path from u to v in Gk, the maximum value that µkv - µku can take subject to constraint (4) in the dual problem, is the shortest path distance from u to v in Gk.

Thus, maximizing µkk - µkr is equivalent to the problem of assigning a length λkuv to each edge in E, subject to constraints (5) and (6) of the dual problem, to maximize the shortest distance from r to k. To maximize the objective of the dual problem, the lengths of {λkuv } are assigned, subject to constraints (5) and (6), to maximize the sum of the shortest distances from r to k for all k ∈ T .

In each iteration (line 5-18), our algorithm increases the dual objective by increasing the value of µkk - µkr for one terminal k that is not connected to r in graph G′ , while maintaining feasibility. It first finds a set of edges δ(k) cutting k and r (line 7). Lines 8-16 calculate a maximum value ∆(s∗ , k) that all λkuv with (u, v) ∈ δ(k) can increase by without violating constraints (5) and (6), which therefore increases the value of µkk - µkr by the same amount. In line 8-10, for each Es, the maximum amount of value that all variables λkuv with (u, v) ∈ Es ∩ δ(k) can increase is equal to the remaining value of constraint (5) divided by the number of edges in Es ∩δ(k). For the edge (r, k) ∈ ET , since it is always in δ(k) and the value of λkr,k can at most increase to 1 by constraint (6), the maximum value ∆(s∗ , k) is at most (1 - λkr,k ) (line 11-12). If ∆(s∗ , k) is equal to (1 - λkr,k ) eventually, the constraint (6) will become tight for (r, k), so zrk is set to be 1 due to the complementary slackness condition 2 and the edge (r, k) is added to G′ (line 13). Otherwise, the constraint (5) will become tight for s∗ , so ys∗ is set to be 1 due to the complementary slackness condition 1 and all edges in Es∗ are added to G′ (lines 15). In line 17, dual variables are updated. The algorithm terminates when all terminals are connected to r. Finally, in line 19, we check whether there exists some purchased Es that can be removed without disconnecting any terminal from r. If we find such Es, we set ys = 0.

4.3 Improvement Techniques

Several techniques can improve the results of Algorithm 1.

Expanding Direction In line 7, the edge sets cutting r and k can be found in two different ways:

• Forward method: The cutset δ(k) contains an edge (u, v) if u is connected to r in the constructed graph G′ and k is connected to v in the original graph G. Namely, the path from r to k is created in a forward direction.
• Backward method: the cutset δ(k) contains an edge (u, v) where k is connected to v in G′ and u is connected to r in G. Namely, the path from r to k is created in a backward direction.

We found that each method may perform better than other depending on the budget, so we use both methods in our implementation.

Local Search In the experiments, we found that different orders of picking terminals can give slightly different solution qualities, so we experimented with a local search method to further improve the solution. We start with the β value returned by the bisection procedure and run the PCSW-DSG solver using a different, randomly selected order. If C(β) < B, we decrease β by a small amount (by Proposition 1). Otherwise, we increase β by a small amount. We repeat this process for 10-20 iterations and in each iteration, we re-select the order of picking terminals randomly. Finally, we return the best feasible solution found. We note, however, that the improvement is limited, indicating that the order of picking terminals doesn't affect the quality that much.

Greedy Padding Since we find a near-optimal β rather than an optimal one, the computed solution may not use all the available budget. Therefore, to improve the solution, a greedy procedure can be used to continue purchasing edge sets until the budget is exhausted. In our experiments, only a few additional edge sets can be purchased so the procedure is very fast, but the improvement is not significant, usually 1-2%.

5 Extensions to Other Problems

BSW-DSG is a very general problem that can capture a range of network optimization problems since the candidate edge sets E can be defined arbitrarily, unlike RCW that is limited to adding patches. For example, our framework is applicable when the underlying network is a road network or river network [Wu et al., 2014a; 2014b] where the investment is used to improve certain edge segments rather than certain locations in the network.

The techniques we propose in this paper can also solve the quota problem (QSW-DSG) by using a slightly different bisection procedure: when narrowing the interval of β, we can check whether Z(β) is greater than the required number of terminals being connected to r and then abandon one half of the interval accordingly. The QSW-DSG problem is an important problem in conservation planning and scheduling. For example, the problem studied by Xue et al. [2012] is a special instance of the QSW-DSG problem where all terminals are required to be connected to r. Xue et al. [2015] mention a useful extension to their problem that only requires a fixed fraction of all terminals to be connected to r, which is equivalent to the QSW-DSG problem. While their paper does not give a solution to that problem, our algorithm can solve it approximately. Furthermore, our algorithm can quickly produce a curve to visualize the tradeoff between cost and the fraction of terminals that need to be connected, by solving the PCSW-DSG problem with different parameters, which offers a useful new capability to end-users.

Optimal solutions to the prize-collecting problem are Pareto optimal with respect to cost and coverage; varying β explores different Pareto optimal solutions. Thus, our algorithm can explore different (approximate) Pareto optimal solutions, but it only partially explores the space and does not guarantee

• (a) Optimality w.r.t. budget

(b) Runtime w.r.t. budget

Figure 4: Optimality and runtime versus budget with 10 samples

generating all Pareto optimal solutions. Finally, our algorithm can be potentially extended to solve the network topology optimization problem in social networks [Khalil et al., 2014]. The applicability of our algorithm to this broad list of problems underscores its significance.

6 Experimental Results

6.1 Empirical Data of RCWProblem

We used the data set for the RCW problem introduced by Sheldon et al. [2010]. The study area includes 2537 patches and 443 parcels. The data set specifies the colonization probabilities between all pairs of patches, the grouping of patches into parcels, the prices of parcels and the initial occupancy status of patches. We used a planning horizon of 100 years and compared our algorithm with a MIP solver and with a greedy algorithm. In each iteration, the greedy algorithm chooses the edge set with the highest ratio of increase in objective value to cost, which gives better solutions than another greedy variant that chooses the edge set with the highest increase in objective value. We used the Gurobi Optimzer as the MIP solver [Gurobi, 2015]. All the experiments were run on a 2.2GHz Intel Core i7 CPU with 16GB of RAM.

Small Sample Size We compare three algorithms for 10 samples and the results are shown in Fig. 4. For our algorithm, we applied Algorithm 1 twice, once using the forward method and once using the backward method, and then applied the local search and greedy padding techniques to the better of those two solutions. In Fig. 4(b), our algorithm runs significantly faster than MIP and 10-20 times faster than the greedy algorithm. In Fig. 4(a), both our algorithm and the greedy algorithm produce near optimal solutions. Our algorithm outperforms the greedy algorithm on larger budget sizes and performs slightly worse on smaller budget sizes.

Large Sample Size We also tested our algorithm and the greedy baseline using 300 samples. The results are show in Fig. 5. We did not test Gurobi because it fails to solve MIPs with more than 20 samples by either using too much time or running out of memory. Our algorithm performs almost identically to the greedy algorithm when the budget is smaller than 3 × 108 and outperforms the greedy algorithm when the budget is greater than 3 × 108 . Fig. 5(b) shows that our algorithm is 10-30 times faster than the greedy algorithm and its running time is almost constant with respect to budget, while the greedy algorithm takes longer for larger budgets.

• (a) Optimality w.r.t. budget

(b) Runtime w.r.t. budget

• (a) Corridor design problem (b) Optimality w.r.t. budget

Figure 6: Experiments on synthetic data

This implies that our algorithm can scale to larger network sizes or more samples, that latter of which leads to improved solution quality within the SAA methodology.

6.2 Synthetic Data

In the RCW problem, even though the discrete optimization problem is NP-hard, the greedy algorithm is still able to produce a solution within 80% of the optimal value, which suggests that the problem does not badly violate submodularity. In this section, we design a group of problems that are more challenging. These synthetic problems are motivated by the corridor design problem, which is another important conservation planning problem [Gomes, 2009] in which the goal is to purchase a subset of parcels to build a (long) corridor connecting distant habitat areas in a fragmented landscape. When formulated within our context, the objective function violates submodularity because purchasing individual parcels has no or little effect on the objective function until the entire corridor is purchased, at which point the endpoints become connected and there is a large jump in the objective value.

With this motivation, we construct a small problem instance as follows (see Fig. 6(a)). We first create M1 parcels, each of which are directly connected to the source population and so carry an immediate reward of p1 if purchased. We also create a free parcel which has reward of p2 if it is reachable from the source. The free parcel is connected by a corridor of M3 parcels, which have reward of p3 if accessible. All parcels in the corridor must be purchased to realize the reward of the free parcel. If we set p2 ≫ p1 > p3 , the optimal strategy should first purchase all parcels in the corridor whenever the budget allows, and first purchase all of the M1 parcels if the budget is not enough to purchase the corridor.

The experimental results of our algorithm and the greedy algorithm on the synthetic data are shown in Fig. 6(b). Our algorithm gives the optimal solution as its curve overlaps with the curve of the MIP solver. When budget is not enough (less than 7500) to purchase all parcels in the corridor, the greedy algorithm performs optimally. But it performs poorly when budget is enough (greater than 7500) because it never purchases any parcels in the corridor before all M1 parcels are purchased. From these results, we observed that when the objective function is not submodular as in our simple problem setting, the greedy algorithm may produce arbitrarily bad solutions, while our algorithm can still produce near-optimal solutions.

7 Conclusion

We propose a fast approximate algorithm to solve a spatial conservation planning problem in which a set of budgetconstrained land parcels are purchased to maximize the expected spread of a population of a species. Our algorithm can approximately solve three different variants of a more general network design problem: the budget, quota, and prizecollecting extensions of the Set-Weighted Directed Steiner Graph Problem. It therefore has potential to apply to a variety of other applications, which we plan to explore in future work. On the Red-cockaded Woodpecker problem, our algorithm produces near optimal solutions and runs significantly faster than a standard mixed integer program solver. Also, the solution quality is comparable to or better than a greedy baseline (which is also near optimal on this problem), but our algorithm is 10-30 times faster. On synthetic problems that do not exhibit submodularity, our algorithm significantly outperforms the greedy baseline.