Optimal Cybersecurity Investments in Large Networks Using SIS Model: Algorithm Design

We study the problem of minimizing the (time) average security costs in large networks/systems comprising many interdependent subsystems, where the state evolution is captured by a susceptible-infected-susceptible (SIS) model. The security costs reflect security investments, economic losses and recovery costs from infections and failures following successful attacks. We show that the resulting optimization problem is nonconvex and propose a suite of algorithms - two based on a convex relaxation, and the other two for finding a local minimizer, based on a reduced gradient method and sequential convex programming. Also, we provide a sufficient condition under which the convex relaxations are exact and, hence, their solution coincides with that of the original problem. Numerical results are provided to validate our analytical results and to demonstrate the effectiveness of the proposed algorithms.


INTRODUCTION
T ODAY, many modern engineered systems, including information and communication networks and power systems, comprise many interdependent systems. For uninterrupted delivery of their services, the comprising systems must work together and oftentimes support each other. Unfortunately, this interdependence among comprising systems also introduces a source of vulnerability in that it is possible for a local failure or infection of a system by malware to spread to other systems, potentially compromising the integrity of the overall system. Analogously, in social networks, contagious diseases often spread from infected individuals to other vulnerable individuals through contacts or physical proximity.
From this viewpoint, it is clear that the underlying networks that govern the interdependence among systems have a large impact on dynamics of the spread of failures or malware infections. Similarly, the topology and contact frequencies among individuals in social networks significantly influence the manner in which diseases spread in societies. Thus, any sound investments in security of complex systems or the control of epidemics should take into account the interdependence in the systems and social contacts in order to maximize potential benefits from the investments.
In our model, attacks targeting the systems arrive according to some (stochastic) process. Successful attacks on the systems can also spread from infected systems to other systems via aforementioned dependence among the systems. The system operator decides appropriate security investments to fend off the attacks, which in turn determine their breach probability, i.e., the probability that they fall victim to attacks and become infected.
Our goal is to minimize the (time) average costs of a system operator managing a large system comprising many Any mention of commercial products in this paper is for information only; it does not imply recommendation or endorsement by NIST. systems, such as large enterprise intranets. The overall costs in our model account for both security investments and recovery/repair costs ensuing infections or failures, which we call infection costs in the paper. To this end, we first consider a scenario where malicious actors launch external or primary attacks. When a primary attack on a system is successful, the infected system can spread it to other systems, which we call secondary attacks by the infected systems, to distinguish them from primary attacks. When primary attacks do not stop, it is in general not possible to achieve an infection-free state at steady state.
In the second case, we assume that there are no primary attacks and examine the steady state, starting with an initial state where some systems are infected. The goal of studying this scenario is to get additional insights into scenarios where the primary attacks occur infrequently. It turns out that, even in the absence of primary attacks, infections may persist due to secondary attacks and the system may not be able to attain the infection-free steady state, and we can compute an upper bound on the optimal value, which is also tight under some condition, more easily.
We formulate the problem of determining the optimal security investments that minimize the average costs as an optimization problem. Unfortunately, this optimization problem is nonconvex and cannot be solved easily. In order to gauge the quality of a feasible solution, we obtain both a lower bound and an upper bound on the optimal value of our problem. A lower bound can be acquired using one of two different convex relaxations of the original problem we propose. For the convex relaxations, we also derive a sufficient condition under which the solution of the convex relaxation solves the original nonconvex optimization problem (Lemma 3). An upper bound on the optimal value can be obtained using an algorithm that finds a local minimizer. Here, we propose two methods -a reduced gradient method (RGM) and sequential convex programming (SCP), both of which produce a local minimizer. Together, our approach offers a bound on the optimality gap.
Numerical studies show that the computational require-ments for the proposed methods are light to modest even for large systems, except for one method, which requires the calculation of an inverse matrix. They suggest that, in almost all cases that we considered, the gap between the lower bound on the optimal value and the cost achieved by our solutions is small; in fact, in most cases, the gap is less than 2-3 percent with the gap being less than 0.3 percent in many cases. In addition, when the infection costs are large, which are likely true in many practical scenarios, the sufficient condition for the convex relaxations to be exact holds, and we obtain optimal points by solving the convex relaxations. Finally, the RGM is computationally most efficient (with the computational time being less than two seconds in all considered cases and less than 0.1 seconds in most cases) and the quality of solutions is on par with that of other methods. This suggests that the RGM may offer a good practical solution for our problem.

Related Literature
Given the importance of cybersecurity, robustness of complex systems, and control of epidemics, there is already a large body of literature that examines how to optimize the (security) investments in complex systems [15], [22], the mitigation of disease or infection spread [7], [11], [30], feasibility and case studies of cyber insurance using pre-screening or differentiated pricing based on the security investments of the insured [16], [33], and designing good attack models and effective mitigating defense against attacks [27], [29], [40]. In view of the volume of existing literature, here we summarize only a small set of studies most closely related to our study.
In [14], [17], [20], [22], the authors adopted a game theoretic formulation to study the problem of security investments with distributed agents or autonomous systems that do not coordinate their efforts. The problem we study in this paper is complementary, but is very different from those studied in the aforementioned studies: in our setting, we assume that the system is managed by a single operator and is interested in minimizing the average (security) costs over time by determining (nearly) optimal security investments. Our study is applicable to, for example, the problem of finding suitable security investments in large enterprise intranets supporting common business processes or supervisory control and data acquisition systems comprising many subsystems.
In another line of research, which is most closely related to our study, researchers investigated optimal strategies using vaccines/immunization (prevention) [7], [36], antidotes or curing rates (recovery) [4], [25], [32] or a combination of both preventive and recovery measures [30], [37]. For example, [36] studies the problem of partial vaccination via investments at each individual to reduce the infection rates, with the aim of maximizing the exponential decay rate to control the spread of an epidemic. Similarly, [25] examines the problem of determining the optimal curing rates for distributed agents under different formulations. In particular, the last formulation of the problem [25], for which only partial result is obtained, is closely related to a special case of our formulation studied in Section 6.
Key differences between existing studies, including those listed above, and ours can be summarized as fol-lows: first, unlike previous studies that focus on either the expected costs from single or cascading failures/infections [19], [20], [22] or the exponential decay rate to the diseasefree state as a key performance metric, we aim to minimize the (time) average costs of a system operator, while accounting for both security investments and infection costs, with both primary and secondary attacks, by modeling time-varying states of systems due to the transmissions of failures/infections. In the presence of primary attacks, it is in general not possible to achieve the infection-free steady state and, thus, the exponential decay rate is no longer a suitable performance metric for our study. Second, unlike some studies that assume that the expected costs/risks seen by systems are convex functions of security investments (e.g., [15]), the expected risks are derived from the steady state equilibrium of a differential system that describes system states and depends on security investments. As it will be clear, the lack of a closed-form expression for the steady state equilibrium complicates considerably the analysis and algorithm design.
Preliminary results of this paper were reported in [26]. In this paper, we extend the findings of [26] in several significant directions. First, we offer an alternative convex relaxation of the original problem. As we demonstrate, this new relaxation technique based on exponential cones, is more efficient and avoids the key issues of the approach based on M-matrix theory presented in [26]. Second, we present another computationally efficient approach to finding a suboptimal solution using SCP together with our Mmatrix theory-based approach, which provides an upper bound on the optimal value of the original nonconvex optimization problem. We compare this approach to one based on the RGM and show that the quality of solutions from these two methods is comparable, but the RGM holds a slight computational edge. Finally, we study a special case with no primary attacks, which is related to the epidemic control problem studied in [25], [28], [32]. Although this can be viewed as a limit case of our problem formulation as the rates of primary attacks go to zero, our approaches to obtaining upper and lower bounds of the optimal value require significant modifications for the reason explained in Section 6. Moreover, we derive sufficient conditions for optimality, which can be verified relatively easily.
The rest of the paper is organized as follows: Section 2 explains the notation and terminology we adopt. Section 3 describes the setup and the problem formulation, including the optimization problem. Section 4 discusses two different convex relaxations of the original problem, followed by two methods for finding local minimizers in Section 5. We discuss a special case with no primary attacks in Section 6. Numerical results are provided in Section 7, followed by a discussion on how our formulation and results can be extended in Section 8. We conclude in Section 9.

Notation and Terminology
Let R and R + denote the set of real numbers and nonnegative real numbers, respectively. Given a set A, we denote the closure, interior, and boundary of A by cl(A), int(A), and ∂A, respectively.
For a matrix A = [a i,j ], let a i,j denote its (i, j) element, A T its transpose, ρ(A) its spectral radius, and σ(A) and σ(A) the smallest and largest real parts of its eigenvalues. For two matrices A and B, we write A ≥ B if A − B is a nonnegative matrix. We use boldface letters to denote vectors, e.g., x = [x 1 , ..., x n ] T and 1 = [1, ..., 1] T . For any two vectors x and y of the same dimension, x • y and x y are their element-wise product and division, respectively. For x ∈ R n , diag(x) ∈ R n×n denotes the diagonal matrix with diagonal elements x 1 , . . . , x n .
A directed graph G = (V, E) consists of a set of nodes V, and a set of directed edges E ⊆ V×V. A directed path is a sequence of edges in the form (i 1 , i 2 ), (i 2 , i 3 ), ..., (i k−1 , i k ) . The graph G is strongly connected if there is a directed path from each node to any other node.

M-Matrix Theory
The set of (nonsingular) n × n M-matrices is denoted by (M n×n + ) M n×n . Note that this definition implies that the off-diagonal elements of A are nonpositive and the diagonal elements are nonnegative; any matrix satisfying these conditions is called a Z-matrix. We shall make use of the following results on the properties of a nonsingular M-matrix [35].
and only if one of the following conditions holds: The next result is a direct consequence of [23,Thm. 2].
i,j is a convex and decreasing function in z ∈ R n + for all 1 ≤ i, j ≤ n.

MODEL AND FORMULATION
Consider a large system consisting of N systems, and denote the set of comprising systems by A := {1, 2, . . . , N }. The security of the systems is interdependent in that the failure or infection of a system can cause that of other systems. 1 As stated before, we study the problem of determining security investments for hardening each system in order to defend the systems against attacks in large systems, in which the comprising systems depend on each other for their function. The goal of the system operator is to minimize the average aggregate costs for all systems (per unit time), which account for both security investments and any economic losses from failures/infections of systems.

Setup
We assume that each system experiences primary attacks from malicious actors. Primary attacks on system i ∈ A occur in accordance with a Poisson process with rate λ i ∈ R + . When a system experiences an attack, it suffers an infection and subsequent economic losses with some probability, called breach probability. This breach probability depends on the security investment on the system: let s i ∈ R + be the security investment on system i (e.g., investments in monitoring and diagnostic tools). The breach probability of system i is determined by some function q i : R + → (0, 1]. In other words, when the operator invests s i on system i, its breach probability is equal to q i (s i ). We assume that q i is decreasing, strictly convex and continuously differentiable for all i ∈ A. It has been shown [2] that, under some conditions, the breach probability is decreasing and log-convex.
When system i falls victim to an attack and becomes infected, the operator incurs costs c r i per unit time for recovery (e.g., inspection and repair of servers). Recovery times are modeled using independent and identically distributed (i.i.d.) exponential random variables with parameter δ i > 0. Besides recovery costs, the infection of system i may cause economic losses if, for example, some servers in system i have to be taken offline for inspection and repair and are inaccessible during the period to other systems that depend on the servers. To model this, we assume that the infection of system i introduces economic losses of c e i per unit time. Besides primary attacks, systems also experience secondary attacks from other infected systems. For example, this can model the spread of virus/malware or failures in complex systems. The rate at which the infection of system i causes that of another system j is denoted by β i,j ∈ R + . When β i,j > 0, we say that system i supports system j or system j depends on system i. Let B = [b i,j : i, j ∈ A] be an N × N matrix that describes the infection rates among systems, where the element b i,j is equal to β j,i . We adopt the convention β i,i = 0 for all i ∈ A.
Define a directed graph G = (A, E), where a directed edge from system i to system j, denoted by (i, j), belongs to the edge set E if and only if β i,j > 0. We assume that matrix B is irreducible. Note that this is equivalent to assuming that the graph G is strongly connected.

Model
We adopt the well-known susceptible-infected-susceptible (SIS) model to capture the evolution of system state. Let p i (t) be the probability that system i is at the 'infected' state (I) at time t ∈ R + . We approximate the dynamics of p(t) := (p i (t) : i ∈ A), t ∈ R + , using the following (Markov) differential equations, which are derived in [28] and are based on mean field approximation. This model is also similar to those employed in [11], [25], [30], [32], [36], [37]: for fixed security investments, s = (s i : i ∈ A) ∈ R N + , In practice, the breach probability q i can be a complicated function of the security investment. Here, in order to make progress, we assume that the breach probability functions can be approximated (in the regime of interest) using a function of the form q i (s) = (1 + κ i s) −1 for all i ∈ A. The parameter κ i > 0 models how quickly the breach probability decreases with security investment for system i. The assumed function satisfies log-convexity shown in [2].
Define α i := κ i δ i , i ∈ A, and α := (α i : i ∈ A). The following theorem tells us that, for a fixed security investment vector s := (s i : i ∈ A) ∈ R N + , there is a unique equilibrium of the differential system described by (1). Theorem 1. Suppose λ 0, δ > 0 and s ≥ 0 are fixed. If the network is strongly connected, i.e., B is irreducible, there exists a unique equilibrium p * ∈ (0, 1) N of (1). Moreover, starting with any p 0 satisfying p * ≤ p 0 ≤ 1, the iteration converges linearly to p * with some rate ρ 0 < 1 − min i∈A p * i . Proof. Please see Appendix A.
Note that the unique equilibrium of the differential system given by (1) specifies the probability that each system will be infected at steady state. For this reason, we take the average cost of the system, denoted by C avg (s), to be , and w(s) quantifies the security investment costs (per unit time), e.g., w(s) = i∈A s i . We assume that w is continuous, (weakly) convex and strictly increasing, and refer to c simply as the infection costs (instead of infection costs per unit time).
A major difficulty in minimizing the average cost in (3) as an objective function is that the equilibrium p * (s) does not have a closed-form expression. As a result, we cannot simply substitute a closed-form expression for the equilibrium p * (s) in (3) and minimize the average cost with s as the optimization variables. For this reason, we formulate the problem of determining optimal security investments that minimize the average cost C avg (s) as follows: where g(s, p) = (g i (s, p) : i ∈ A), and Recall that, for given s ∈ R N + , only the unique equilibrium p * ∈ (0, 1) N in Theorem 1 satisfies the constraint in (4b).
Clearly, the solution to problem (P) will also shed light on which systems are more critical from the security perspective and, hence, should be protected. In the problem (P), we do not explicitly model any total budget constraint on security investments for simplicity of exposition. However, we will revisit the issue of constraints on security investments, such as a total budget constraint, and discuss how it affects our main results in Section 8.1.
This problem (P) is nonconvex due to the nonconvexity of the equality constraint functions in (4b). In particular, g i contains both quadratic or bilinear terms p i p j and p i s i . In the following sections, we develop four complementary algorithms for finding good-quality solutions to the nonconvex problem: the first two approaches are based on a convex relaxation using different techniques, and provide a lower bound on the optimal value of the problem (P). The last two are designed to find a local minimizer of the problem (P), hence provide an upper bound on the optimal value, and are based on the RGM and SCP.

LOWER BOUNDS VIA CONVEX RELAXATIONS
In this section, we discuss how we can relax the original problem (P) and construct two different convex formulations, which can be used to obtain (a) a lower bound on the optimal value and (b) a feasible solution to (P) using optimal points of the relaxed problems. Furthermore, we provide a sufficient condition for the relaxed problems to be exact, i.e., their optimal point is also an optimal point of the nonconvex problem (P). The first approach is based on M-matrix theory and the preliminary results were reported in [26]. The second approach is designed to deal with some of computational issues of the first approach. Moreover, as we will show, the optimal point of the first approach can be computed from that of the second approach.

Convex Relaxation: M-Matrix Theory
Given λ 0 and irreducible B, Theorem 1 states that the unique equilibrium of (1) which satisfies (4b) is strictly positive. Hence, we can rewrite the constraints in (4b) as where p −1 = (p −1 i : i ∈ A). By introducing a new variable the constraint in (5) can be rewritten as Note that (7) is affine in z, s and p, and the nonconvexity in the equality constraint functions (mentioned at the end of the previous section) is now captured by z, which from (6) can be expressed as We can show that the matrix (diag(z)−B) is a nonsingular M-matrix and, hence, p = (diag(z) − B) −1 λ as follows: from (8), since λ 0, we have λ i * > 0 for some i * . Since matrix B is assumed irreducible, for any j such that λ j = 0, we can find a finite sequence Because p > 0, Lemma 1-(f) tells us that this is equivalent to matrix (diag(z) − B) being a nonsingular M-matrix. As a result, the original problem (P) can be reformulated as where We can show that the set Ω in (9) is convex. This is proved in Appendix B. Also, it follows from Lemma 2 that for any 1 ≤ i, j ≤ N , the element diag(z) − B −1 i,j is convex and (element-wise) decreasing in z ∈ Ω. For these reasons, we obtain the following convex relaxation of (P2).
This convex relaxation can be solved by numerical convex solvers to provide a lower bound on the optimal value of (P). Also, as shown in the following theorem, its optimal point also leads to a feasible point for problem (P).
denote an optimal point of (P R1 ) and f * the optimal value of (P). Then, we have . Also, it satisfies (10a) with equality. Thus, it is a feasible point for problem (P), proving the second inequality.
). Based on this, we can provide a following sufficient condition for convex relaxation (P R1 ) to be exact. Lemma 3. The above convex relaxation (P R1 ) is exact if Proof. Suppose (s,p) is the feasible point of (P) given in (11), together with p * R ≥p, this gap is nonpositive. By Theorem 2, this gap must be zero. Thus, (P R1 ) is exact. Remark 1. (Sufficient condition for exact relaxation) First, roughly speaking, the condition in (11) means that when the infection costs c are sufficiently high, the convex relaxation (P R1 ) is exact and we can find optimal security investments, i.e., a solution to (P), by solving (P R1 ) instead. The intuition behind this observation is the following: as c becomes larger, the second term in the objective function, namely c T p, becomes more important and an optimal point tries to suppress it by reducing p. However, since p must satisfy the inequality in (10a), it can only be reduced till the equality holds, which satisfies the constraint in problem (P2). Second, condition (11) can be verified prior to solving the relaxed problem. This can be done easily if w is a linear function or an upper bound on the gradient ∇w is known. Finally, even when the convex relaxation is not exact, (s,p) can still be used as a good initial point for a local search algorithm, such as the RGM developed in Section 5 below.
Remark 2. (Numerical issues of (P R1 )) Although (P R1 ) is a convex problem, there are a few numerical challenges. First, the Jacobian of constraint functions in (10a), which involves the derivative of inverse matrix (diag(z) − B) −1 , tends to be dense even when B is sparse. Thus, off-the-shelf convex solvers may not be suitable for large systems.
Second, although the constraint set Ω for z (defined in (9)) is convex, it is not numerically easy to handle, especially for large networks. This is because Ω is not closed and (P R1 ) becomes invalid outside Ω. Thus, a numerical algorithm ought to stay inside Ω and, for this reason, the nonsingularity of the M-matrix, diag(z) − B, should be ensured at every step. In general, it takes O(N 3 ) to check if the matrix satisfies this condition [34]. The following approach can, however, alleviate the computational burden.
, the solution is optimal for (P R1 ); otherwise, construct a newΩ so that z * R belongs to the interior of newΩ and repeat. Below, we propose an efficient way to choose the subsetΩ that is more suitable for numerical algorithms.

Construction of Convex Subsets of Ω
A key observation to constructing suitable subsets of Ω is that, in view of Lemmas 1 and 2, Ω can be expressed as Thus, for every z ∈ ∂Ω,Ω(z) := {z ∈ R N + | z z} ⊂ Ω. Our goal is to find somež ∈ ∂Ω such that an optimal point x * R that solves the relaxed problem with Ω replaced byΩ(ž), satisfies z * R ∈ intΩ(ž). Below, we provide several possible choices for z with increasing computational complexity.
This can be guaranteed by choosing z > B1, where the lower bound B1 represents the total rate of infection from immediate neighbors in the graph G. From (10b), a trivial sufficient condition is δ + λ ≥ B1. But, we observe empirically that this often leads to suboptimal solutions. 4.1.1.2 Dominant eigenvalue: Another straightforward lower bound is given by z > ρ(B)1. Recall that the spectral radius ρ(B) is also an eigenvalue of B and equal tō σ(B), which can be computed efficiently using, for example, the power method.
4.1.1.3 Iterative dominant eigenvalue selection via matrix balancing: Unfortunately, we observe empirically that a static selection of the subsetΩ does not always lead to a good solution and a following iterative algorithm yields better performance: let h > 0 be a normal vector of the plane tangent to the closure of Ω at some z ∈ ∂Ω such that (12) where the second equality follows from the fact that we are minimizing a linear function over a closed convex set. The minimization in (12) amounts to finding the smallest diagonal perturbation z (in 1-norm weighted by h) so that B becomes (negative) stable. In Appendix C, we show that this is in fact a matrix balancing problem, for which efficient algorithms exist (see [5], [31] for nearly-linear time centralized algorithms and [24], [25] for distributed algorithms with geometric convergence).
Our first proposed algorithm (Algorithm 1) is based on the discussion in this subsection. Initially, we choose somē h > 1 and h = α −1 • ∇w(s 0 ), where s 0 is the initial choice of security investments. This heuristic is based on the relaxed problem by weighting only the investment cost w(s) without considering p. Subsequent iterations are based on dominant eigenvalues with varying weights determined by h + , which reflects active constraints ofz R (of the current solution). Sincez R ∈ Ω, we have σ diag(z R ) − B > 0. Thus, we can construct a new subsetΩ(z) by translating the set {z ≥z R } towards the boundary ∂Ω in the direction of h + , so thatz R lies in the interior of newΩ(z). In our numerical studies (Section 7), we useh = 10.
Note that Algorithm 1 is guaranteed to converge because the problem is convex and the objective function value decreases after each iteration. Although we cannot provide a convergence rate, numerical studies in Section 7 show that only a few iterations are needed in most cases.

Convex Relaxation Based on Exponential Cones
As explained in the previous subsection, a possible difficulty in solving the convex relaxation in (P R1 ) is taking into account two constraints -constraint in (10a) and z ∈ Ω. Here, we present an alternative convex relaxation of the original problem, which avoids these issues by introducing auxiliary optimization variables and relaxing the equality constraint in (4b) without the need for the constraint set Ω.
First, recall from the previous subsection that the constraint in (4b) can be rewritten as Since any solution must satisfy p ∈ (0, 1] N , we introduce following auxiliary variables and rewrite the equality constraint in (5): for fixed y ∈ R N + , define p := e −y , t := λ • e y , U := diag(e y )Bdiag(e −y ). (13) Using these new variables, (5) can be rewritten as follows.
Then, problem (P) is equivalent to the following problem.
The equivalent problem (P3) is still nonconvex due to the constraints in (13). We can relax these equality constraints with the following inequality convex constraints.
This leads to the following second convex relaxation. (14), (15) We can express the constraints in (15) as a following set of at most 2N + m exponential cone constraints: These constraints can be handled efficiently by conic optimization solvers, e.g., MOSEK [1].
Remark 3. We demonstrate below that, somewhat surprisingly, the convex relaxations in (P R1 ) and (P R2 ) are in fact equivalent. Moreover, although one may suspect that the size of (P R2 ) with 4N + m variables and 3N + m constraints is much larger than the size of (P R1 ), the constraints of (P R2 ) are much easier to handle numerically. We will provide numerical results to illustrate this in Section 7.
Analogously to Theorem 2, the following theorem tells us how to find a feasible point of the problem (P), using an optimal point of problem (P R2 ). In addition, it asserts that the two convex relaxations (P R1 ) and (P R2 ) are equivalent in that their optimal values coincide and we can find an optimal point of (P R1 ) from an optimal point of (P R2 ) .
Proof. Please see Appendix D.
As a direct consequence of the theorem, the suboptimality of (s ′ , p ′ ) can be assessed using the gap f (s ′ , p ′ ) − f (s + , p + ). Similar to Lemma 3, the condition in (11) provides a sufficient condition for this gap to be zero, i.e., the convex relaxation in (P R2 ) is exact.

UPPER BOUNDS ON OPTIMAL VALUE
The previous section described (i) how we can formulate a convex relaxation of problem (P), which provides a lower bound on the optimal value of (P), using two different techniques and (ii) how to find a feasible solution to (P) using an optimal point of a convex relaxation.
Although the convex relaxation (P R1 ) or (P R2 ) may be exact under certain conditions, this is not true in general. In addition, (P R1 ) may not scale well due to the constraint in (10a); see also Remark 2 above and numerical results in Section 7. For these reasons, we also propose efficient algorithms for finding a local minimizer of the nonconvex problem (P) in this section. These algorithms provide an upper bound on the optimal value, which, together with the optimal value of a convex relaxation when available, can be used to offer a bound on the optimality gap.

Reduced Gradient Method
Among different nonconvex optimization approaches, we first choose the RGM [10], [21] because it is well suited to the problem (P) and, more importantly, is scalable.

Main Algorithm
First, together with Theorem 1, the implicit function theorem tells us that the condition g(s, p) = 0 in (4b) defines a continuous mapping p * : s ∈ R N + → p * (s) ∈ (0, 1) N such that g(s, p * (s)) = 0. Thus, problem (P) can be transformed to a reduced problem only with optimization variables s: Suppose that (s ⋆ , p ⋆ ) is a feasible point of (P). Then, the gradient of F at s ⋆ is equal to This matrix can be computed by totally differentiating g(s, p * (s)) = 0 at s ⋆ : the calculation of total derivative yields The following lemma shows that M (s ⋆ ) is nonsingular. (20) and the gradient of F is given by Hence, we can apply the (projected) gradient descent method on the reduced problem in (19). For instance, [3,Proposition 2.3.3] shows that this method converges to a stationary point under step sizes {γ t } t≥0 chosen by the Armijo backtracking line search. Note that, after each update of s during a search, we need to compute the corresponding p so that (s, p) is feasible for the problem (P). As mentioned earlier, this can be done by using the fixed point iteration in (2).
Our proposed algorithm is provided in Algorithm 2.
Algorithm 2: Reduced Gradient Method

Computational Complexity and Issues
For large systems, a naive evaluation of the gradient ∇F , which requires the inverse matrix M (s ⋆ ) −T , becomes computationally expensive, if not infeasible. For this reason, we develop an efficient subroutine for computing ∇F . This is possible because our algorithm only requires u (in line 4 For fixed t ∈ N := {0, 1, . . .}, the vector u is a solution to a set of linear equations M T u = c, where the matrix M T tends to be sparse for most real graphs G. Thus, there are several efficient algorithms for solving them. In this paper, we employ the power method: let M = D − E, where D and E denote the diagonal part and off-diagonal part of M , respectively. Then, the linear equations are equivalent to c = Du − E T u. Since D is invertible, the following fixed point relation holds: As Hence, the iteration u k+1 = G(u k ) converges to the solution u exponentially fast. Moreover, this iteration is highly scalable because E = diag(1−p)B is sparse, requiring only O(|E|) memory space and O(|E|) operations per iteration.

Sequential Convex Programming Method
In this subsection, we will develop a second efficient algorithm for finding a local minimizer of the original nonconvex problem (P). This will provide another upper bound on the optimal value we can use to provide an optimality gap together with a lower bound from convex relaxations.
Our algorithm is based on SCP applied to the original formulation in (4a)-(4b). To this end, we successively convexify the constraint (4b) using first order approximations. The novelty of our approach is to linearize (only) the terms p i p j and then employ either the convexity result in Lemma 2 or the exponential cone formulation as in Section 4.
At each iteration t ∈ N, we replace the terms p i p j with their first order Taylor expansion at p (t) , resulting in the following partially linearized equality constraint functions: The partial linearization error is equal to g (t) (s, p) − g(s, p) = (p−p (t) )•B(p−p (t) ). When p is close to p (t) , this error will likely be 'small' and we expect the linearization step to be acceptable. This allows us to approximate the constraint (4b) with g (t) (s, p) = 0, which can be rewritten as where λ (t) = λ + p (t) • (Bp (t) ), and This gives us the following subproblem we need to solve at each iteration t ∈ N.
(S1) min s≥0,p≥0 The solution at the t-th iteration is then used to construct a new constraint for the (t + 1)-th iteration, and we repeat this procedure until some stopping condition is met. Unfortunately, the problem in (24) is still nonconvex. But, as we show below, under certain conditions, it can be transformed to a convex problem, which can be solved efficiently.

Convex Formulation Based on M-matrix
First, we show that if (s, p) is feasible for the subproblem in (24), then L (t) +diag(α•s) ∈ M N ×N + : note that this is always a Z-matrix. Also, from (22) and (23), any feasible p for the subproblem must be positive, and given p (t) > 0, we have λ (t) > 0 from its definition. Because p and λ (t) are positive, together with condition (23), Lemma 1-(g) implies that L (t) + diag(α • s) ∈ M N ×N + . This in turn tells us from Lemma 1-(b) that its inverse exists and is nonnegative. Consequently, This allows us to reformulate the subproblem in (24) as follows: we replace p with the above expression in (25) and introduce a new constraint that s belongs to a feasible set Note that Ω (t) is convex (the proof follows similar arguments in Appendix B). The partially linearized subproblem in (24) can now be written as where ζ(s) := c T L (t) + diag(α • s) −1 λ (t) is a convex and decreasing function on Ω (t) in view of Lemma 2. Thus, the problem (P L ) is convex at every iteration t. Note that, when solving (P L ), we need to ensure the constraint s ∈ Ω (t) is satisfied. This can be done in a manner similar to that discussed in subsection 4.1.1.
After computing an optimal point of (P L ) at the t-th iteration, which we denote by s (t+1) , we then find p (t+1) satisfying g(s (t+1) , p (t+1) ) = 0 (constraint (4b)) using the fixed point iteration in (2). Thus, we obtain a feasible solution (s (t+1) , p (t+1) ) to the original problem (P) after each iteration t ∈ N. The proposed algorithm based on this approach is provided in Algorithm 3 below.

Algorithm 3: Sequential Convex Programming
Remark 4. (Complexity) For small and medium-sized networks, the subproblem can be solved using off-the-shelf numerical convex solvers, e.g., interior point methods. In this paper, we use an interior-point method to solve the subproblem (P L ), which employs the Newton's algorithm on a sequence of equality constrained problems. Since the number of variables is O(N ) and the number of constraints is also O(N ), the worst case complexity is O(N 3 ) [5].
For large networks, we take advantage of the fact that we do not need to solve (P L ) exactly at each iteration. Thus, we can use simple approximations of Ω (t) and employ computationally cheaper methods to solve (P L ). For example, we can follow the same gradient-based approach in [23] for solving the subproblem, where the gradient of ζ(s) given by where S = L (t) + diag(α • s), can be computed efficiently using the power method as explained in subsection 5.1.2.

Convex Formulation Based on Exponential Cones
As discussed above, if (s, p) is feasible for problem (S1), L (t) + diag(α • s) is a nonsingular M-matrix and p > 0; see also (25). Thus, we can introduce a new variable y satisfying Then, (23) becomes L (t) +diag(α • s) e −y = λ (t) , which, after left-multiplying both sides by diag(e y ), is equivalent to where B (t) = diag(1 − p (t) )B. As a result, the subproblem is equivalent to the following problem.
(S2) min s≥0,y,t,U A convex relaxation of (S2) can be obtained by replacing the last two equality constraints with inequality constraints.
It turns out that this convex relaxation is always exact. Proof. A proof can be found in Appendix E.
We end this subsection by noting that this convex formulation is in fact equivalent to the one based on M-matrix in (26). This can be shown using similar arguments used in the proof of Theorem 3 and is omitted here.

SPECIAL CASE: λ = 0
In practice, we expect that the systems experience primary attacks infrequently and λ is small, and that steady-state infection probabilities are not large. For this reason, we consider a limit case of our problem as λ → 0 with diminishing primary attack rates. As we show, studying the special case with λ = 0 reveals additional insights into the steady-state behavior and provides an approximate upper bound on the system cost when λ ≈ 0, which can be computed easily.
This case reduces to a problem that has been studied by previous works, in which the adjustable curing rate is equal to δ i /q i (s i ) for each i ∈ A. 2 A key difference between this case and when λ 0 is that Theorem 1 cannot be applied to guarantee the uniqueness of an equilibrium because the assumption λ 0 is violated. It turns out that this difference has significant effects on our problem, as it will be clear.

Preliminary
In the absence of primary attacks, if no system is infected at the beginning, obviously they will remain at the state. However, if some systems are infected initially, there are two possible outcomes based on the security investments s.
Case 2: ρ diag(α • s + δ) −1 B > 1 -In this case, there are two equilibria of (1) -one stable equilibrium p se (s) > 0 and one unstable equilibrium 0: (a) if p(0) = 0, although there are no primary attacks, we have p(t) → p se (s). As a result, somewhat surprisingly, infections continue to transmit among the systems indefinitely and do not go away; and (b) if p(0) = 0, obviously p(t) = 0 for all t ∈ R + .
Based on this observation, we define a function p se : where p se (s) is the aforementioned stable equilibrium of (1) for the given security investment vector s ∈ R N + . It is shown in Appendix I that p se is a continuous 2. In the previous studies [25], [28], [32], security investments affect the curing rates rather than the breach probability, i.e., they determine how quickly each system can recover from an infection, but do not change the infection probability of systems. function over R N + . This tells us that, if we start with p(0) = 0, for any given security investments s ≥ 0, our steadystate cost is given by w(s) + c T p se (s). For this reason, we are interested in the following optimization problem. min s≥0 w(s) + c T p se (s) (28) We denote the optimal value and the optimal set of (28) by f * 0 and S * 0 , respectively. Based on the discussion, we have the following simple observation.
The theorem follows directly from the above discussion and the monotonicity of the spectral radius of nonnegative matrices [13,Thm 8.1.18] Remark 5. The theorem rules out the case where ρ diag(α • s * + δ) −1 B < 1 for some s * = 0. It suggests that if the recovery rates of all the systems are sufficiently large, no additional investments are needed. Otherwise, at any solution s * ∈ S * 0 , the spectral radius is either (i) at the threshold (of one) or (ii) strictly above the threshold. In case (i), w(s * ) is also the smallest investment cost to suppress the spread in that lim t→∞ p(t) = 0 for all p(0). We will show in subsection 6.2 below how to compute this minimum investment for suppression, denoted by C * . On the other hand, case (ii) corresponds to an endemic state, i.e., lim t→∞ p(t) = p se (s * ) > 0 when p(0) = 0. In this case, we will demonstrate that we can find upper and lower bounds on the optimal cost using our techniques in Sections 4 and 5.
In order to facilitate our discussion, we introduce following related optimization problems with a fictitious constraint on security investments. The goal of imposing a fictitious budget constraint is not to investigate a problem with a budget constraint; instead, it is used to facilitate the determination of a (nearly) optimal point of (28) as we will show.
We define a function f 0 : R + → R + , where f 0 (C) is the optimal value of the above optimization problem for a given budget C. Clearly, the function f 0 is continuous and nonincreasing, and lim C→∞ f 0 (C) = f * 0 , i.e., problem (29) reduces to (28) by letting C → ∞.
Obviously, w * is the minimum security investments necessary to minimize the total cost in (28), and s * is an optimal point with the smallest security investments. Then, because f 0 is nonincreasing, for any C ≥ w * , we have which implies f 0 (C) = f 0 (w * ) = f * 0 . On the other hand, where the strict inequality follows from the definition of s * in (30); any s with w(s) < w * is not an optimal point and, as a result, we have f * 0 < w(s) + c T p se (s).
The inequalities in (31) and (32) tell us the following: increasing the security budget C reduces the total cost f 0 (C) while C ≤ w * . On the other hand, beyond w * , increasing the budget will not reduce the cost any more as f 0 (w * ) = f * 0 . (28) From the discussion at the beginning of subsection 6.1, it is clear that the spectral radius of the matrix diag(α • s + δ) −1 B plays an important role in the dynamics and the determination of a stable equilibrium of (1). For this reason, we find it convenient to define the following problem:

Bounds on the Optimal Value of
Let C * be the optimal value of this optimization problem, which is the aforementioned minimum investments needed for suppression. We show in Appendix H that this optimization problem can be transformed into a convex (exponential cone) problem and, thus, can be solved efficiently.
The following lemma points out an important fact that we will make use of in the remainder of the section.
Lemma 5. The optimal value C * of (33) is an upper bound on f * 0 , i.e., f * 0 ≤ C * . Proof. We know that any optimal points 0 of (33) satisfies w(s 0 ) = C * and p se (s 0 ) = 0. Therefore, the total cost achieved bys 0 is equal to From this lemma and the definition of s * in (30), we have Obviously, this also implies f * 0 = f 0 (C) for all C ≥ C * . In a special case when the equalities in (34) hold, we have f * 0 = w * = C * , and an optimal point of the optimization problem in (33), says 0 , is optimal for the problem in (28) with p se (s 0 ) = 0. Also, ρ diag(α • s * + δ) −1 B = 1 and the equality holds in the second part of Theorem 5. But, in general these equalities may not hold, in which case we must have w * < f * 0 < C * and ρ diag(α • s * + δ) −1 B > 1. Because w * is unknown beforehand, we cannot determine which case holds. However, if we solve the constrained problem in (29) with C = C * − ǫ for small positive ǫ, either (a) the optimal point s * C we obtain is an optimal point of (28) if w(s * Note that in the latter case, we have C * −f * 0 ≤ ǫ, ands 0 is ǫ-suboptimal for (28).
At first glance, one may suspect that solving (29) is as difficult as solving (28) because both share the same nonconvex objective function. However, (29) enjoys a few numerical advantages: first, unlike (28), the constrained problem (29) with w(s) ≤ C < C * admits only p se (s) > 0, which can be computed using the fixed point iteration in (2). Second, perhaps more importantly, because the stable equilibrium is guaranteed to be strictly positive, it allows us to employ the approaches in Sections 4 and 5 in order to compute upper and lower bounds on the optimal value f 0 (C). In the process of finding these bounds, we also demonstrate an interesting observation that we can bound the gap C * − f * 0 under certain conditions. Before we proceed with discussing the bounds, let us remark on the choice of ǫ. Theoretically, we want ǫ as small as possible because it determines the ǫ-suboptimality ofs 0 when the constraint is active (case (b) above). In practice, however, ǫ should not be too small because it can cause numerical issues when C * = w * and p se (s * ) = 0 (or p se (s * ) ≈ 0 when C * ≈ w * ). This is because our approaches to finding upper and lower bounds on the optimal cost in Sections 4 and 5 rely on the strict condition that p se (s) > 0.

Lower Bound via Convex Relaxation
Borrowing a similar approach used in subsection 4.2, we can formulate the following convex relaxation of (29): We denote the optimal value of (P R3 ) by f L (C). The inequalities in (35b) and (35c) can be expressed as exponential cone constraints as done in (16a) and (16c), respectively. Thus, this convex relaxation can be solved efficiently.

Together with the earlier inequality in
More can be said regarding these bounds as follows.
Theorem 6. Suppose C = C * − ǫ and (s L , p L , y L , U L ) is an optimal point of (P R3 ). Let C L = w(s L ). Then, we have Proof. Please see Appendix F.
This result tells us that, when C L = C, any optimal point of (33) is O(ǫ)-suboptimal for the original problem in (28). Also, since f L (0) ≤ f 0 (0) = c T p se (0), the bound in (37) is upper bounded by ǫ 1 + c T pse(0) C * , which can be computed before solving (P R3 ). Therefore, a natural question that arises is: Can we determine if the condition C L = C holds for some C < C * without having to solve the convex relaxation (P R3 )?
The following theorem offers a (partial) answer to this question by providing a sufficient condition for the condition C L = C to hold. For a given budget constraint C ∈ R + , define S C = {s ∈ R N + | w(s) ≤ C}. Theorem 7. Suppose that every s ∈ S C satisfies Then, C L = C. If (38) holds for all s ∈ S C * , then f * 0 = C * .
Proof. A proof can be found in Appendix G.
Note that the condition (38) can be verified prior to solving the relaxed problem (P R3 ). In the case that condition (38) holds for all s ∈ S C * , any optimal point s 0 of (33) is also optimal for our original problem in (28).

Upper Bound via the Reduced Gradient Method
In order to use the RGM for the problem in (29), we first need to introduce a following modification to Algorithm 2: replace line 5 of Algorithm 2 with where P SC [·] denotes the Euclidean projection onto S C . When w is simple, this projection step can be very efficient. The results in Section 5 still hold in this case. In particular, at any feasible point (s ⋆ , p ⋆ ) of problem (P) with λ = 0 such that s ⋆ ∈ S C , the matrix which arises from totally differentiating the constraint g(s, p) = 0, is still a nonsingular M-matrix. This can be verified by noting that M (s ⋆ ) satisfies M (s ⋆ )p ⋆ = p ⋆ • (Bp ⋆ ) > 0, where the positivity follows from p ⋆ > 0 because we require that s ⋆ ∈ S C with C < C * . As a result, we can use Algorithm 2 with an efficient evaluation of reduced gradient as shown in subsection 5.1.2 and the projection step as described above.
Remark 6. We summarize how to find a good solution to (28) when ρ diag(δ) −1 B > 1 based on the above discussion: first, find a pair (s 0 , C * ) of the optimal point and optimal value of (33). If (38) holds for all s ∈ S C * , s 0 is an optimal point of (28). Otherwise, solve (P R3 ) with C = C * − ǫ (for small ǫ) and let (s L , f L ) be the pair of its optimal point and optimal value. If w(s L ) = C, adopt s 0 as a solution to (28) with opt_gap ≤ ǫ(1+ c T pse(0) C * ). Otherwise, solve (29) using RGM and adopt its solution s U as a solution to (28) with

NUMERICAL RESULTS
In this section, we provide some numerical results that demonstrate the performance of the proposed algorithms. Our numerical studies are carried out in MATLAB (version 9.5) on a laptop with 8GB RAM and a 2.4GHz Intel Core i5 processor. We consider 5 different strongly connected scale-free networks with the power law parameter for node degrees set to 1.5, and the minimum and maximum node degrees equal to 2 and ⌈3 log N ⌉, respectively, in order to ensure network connectivity with high probability.
For all considered networks, we fix α i = 1 and δ i = 0.1 for all i ∈ A. The infection rates β j,i , (j, i) ∈ E, are modeled using i.i.d. Uniform(0,1) random variables. We choose w(s) = 1 T s and c = νB T 1 + 2c rand , where the elements of c rand are given by i.i.d. Uniform(0,1) random variables and ν ≥ 0 is a varying parameter. We select c above, in order to reflect an observation that nodes which support more neighbors should, on the average, have larger economic costs modeled by c e i (Section 3-A). We consider two separate cases: λ > 0 and λ = 0.

Case λ > 0
We generate λ using i.i.d. Uniform(0,1) random variables for each network, set ν ∈ {0, 0.5, 1}, and apply 5 schemes described below. The results (averaged over 10 runs) are summarized in Table 1, and a more detailed description of the simulation setups can be found in Appendix J. Here, the reported optimality gap is the relative optimality gap given by opt gap are the optimal values of (P R1 ) and (P R2 ), respectively, and f cur is the cost achieved by the solution found by the algorithm under consideration. Although the two convex relaxations are equivalent, their numerical solutions are not necessarily identical and we take a conservative lower bound given by the minimum of the two values. When the optimal value of a convex relaxation is unavailable, we take the other optimal value. Also, the column t s indicates the total runtime.
• M-matrix + OPTI: We solve the relaxed problem based on M-matrix in subsection 4.1 using Algorithm 1, where line 3 utilizes an interior point method from the OPTI package [8], and consider the feasible point (s,p) in Theorem 2. The column iter shows the pair of (i) the number of outer updates (each corresponding to an approximatioñ Ω(z (t) ) of the set Ω) and (ii) the average number of inner interior-point iterations inside outer updates. As we can see, the algorithm runtime does not scale well with the network size; for the case (N, |E|) = (2001, 12076), the solver failed to converge within an hour.
• K-Exp + MOSEK: We solve the relaxed problem (P R2 ) with exponential cone constraints using the MOSEK package [1] and consider the feasible point (s ′ , p ′ ) in Theorem 3. The column iter indicates the number of interior point iterations. As expected, this method enjoys smaller runtimes and, hence, has a computational advantage over the M-matrix + OPTI scheme.
• RGM + ARMIJO: We use the RGM in Algorithm 2 to find a local minimizer (s ⋆ , p ⋆ ). The column iter shows the pair of (a) the number of gradient updates and (b) the maximum number of fixed point iterations needed for evaluating u in line 4 and p * in line 7, denoted byk fp . In our study, the reported values ofk fp are all relatively small as expected from our earlier discussions (Theorem 1 and subsection 5.1.2).
• M-Matrix SCP: We use Algorithm 3 to find a local minimizer. Here, we solve the convex optimization subproblem in line 5 of Algorithm 3 approximately, using the OPTI package for N ≤ 10 3 and a gradient descent method for N > 10 3 . The column iter shows the pair of (a) the number of outer linearization updates and (b) the average inner steps of either the interior-point solver or gradient descent method. However, this approach does not scale well due to the M-matrix based relaxation as shown in Table 1.
• K-Exp SCP: For this algorithm, we replace the convex optimization subproblem in line 5 of Algorithm 3 with the formulation in (27) and solve it approximately using MOSEK. The column iter shows the pair of (a) the number of outer linearization updates and (b) the average inner interior-point steps. As we can see from Table 1, this approach achieves similar opt gap as RGM + ARMIJO, but the runtime is roughly two orders higher. Compared to M-Matrix SCP, its performance, both in terms of the quality of solution and runtime, is comparable.
We summarize observations. First, as ν increases and infection costs become larger, as expected from Lemma 3, the gap diminishes and becomes negligible when ν = 1.  Second, the upper bounds from local minimizers are very close to the lower bound f * R , even when the relaxation may not be exact (for ν = 0 and 0.5). Moreover, they lead to optimal points when the relaxation is exact. This suggests that the algorithms can practically find global solutions to the original problem. Finally, Algorithm 2 based on RGM, is highly scalable: despite a larger number of required iterations compared to all other schemes, the total runtime t s is much smaller and is a fraction of that of Algorithm 1 or 3; we note that we did not optimize step sizes; we instead used the same parameters in all cases.
We also tried sqp and interior-point solvers in MATLAB for problem (P), but found them to be very inefficient compared to our approaches to finding local optimizers. For example, for the case (N, |E|) = (999, 5750) and ν = 0.5, while RGM runs in only a fraction of a second, sqp takes 19 iterations in 102 seconds to achieve the same opt gap as RGM, and interior-point terminates after 125 iterations in 68 seconds with twice the opt gap of RGM.

Case λ = 0
In this subsection, we study the scenario with no primary attacks, using the scale-free network with 499 nodes from the previous subsection. We consider the value of ν in {0.6, 0.8, 1} in order to obtain more informative numbers.
Following the steps outlined in Section 6, we first find the optimal value C * and an optimal point s 0 of problem (33), using MOSEK. When ν = 1, condition (38) in Theorem 7 holds for all s ∈ R N + and, thus, we have f * 0 = C * with s 0 being an optimal point of the original problem in (28).
Second, for ν < 1, we consider the problem in (29) with ǫ = 0.01C * or, equivalently, C 1 = 0.99C * , and find a lower bound f L (C 1 ), which is the optimal value of the relaxed problem (P R3 ), using MOSEK. For ν = 0.6 and 0.8, we found that the constraint w(s) ≤ C 1 is inactive at the optimal point and, consequently, f * L = f L (C 1 ). In addition, using the projected RGM described in subsection 6.2.2, we also compute an upper bound f U (C 1 ) on the optimal value f * 0 and then consider the gap ∆f 0 := min{f U (C 1 ), C * }−f * L . We plot in Fig. 1 both the upper bound f U (C) and the lower bound f L (C) of the optimal value of problem in (29) as a function of C over the interval [0, 0.99C * ]. There are several observations we make from the plots. o1. When ν = 0.6 and the infection costs c are small, the left plot shows f U (C) < C * , which tells us that the local minimizer found by the projected RGM of the problem in (29) with the given security budget C, is better than the optimal point of (33). The plot also indicates that both the upper and lower bounds quickly reach a plateau less than C * with increasing C. This likely suggests that the optimal security investment at an optimal point of (28) is significantly smaller than C * . o2. As we discussed just before subsection 6.2.1, the left plot highlights the practical usefulness of our approach: by considering the problem in (29) for C = C * − ǫ with small positive ǫ, we can quickly estimate the optimal security investments, i.e., whether or not they are close to C * (when the security budget constraint is active at s * C ) or is equal to w(s * C ) (when the constraint is inactive at s * C ), without suffering from the numerical issues mentioned earlier. o3. As ν increases, so do the bounds f U (C 1 ) and f L (C 1 ) (normalized by C * ). For larger values of ν with higher infection costs (middle and right plots), the upper bound f U (C 1 ) obtained from a local minimizer is slightly larger than another upper bound C * , suggesting that the budget constraint is active at a local minimizer returned by the projected RGM. This suggests that, when infection costs are high, s 0 , which may overinvest compared to an optimal point, may still be a good feasible point. Furthermore, as ν gets close to one, eventually s 0 becomes an optimal point for the problem in (28) as shown in the right plot for ν = 1. This is expected because as the infection costs become larger, the system operator has an incentive to invest more in security. o4. Although we do not report detailed numbers here, both the upper and lower bounds can be computed efficiently using RGM and MOSEK; the required computational time is always less than 2 seconds for each run, suggesting that the RGM may be a good method for identifying suitable security investments for large systems.

Constraints on Security Investments
As mentioned in Section 3, we imposed only non-negativity constraints on security investments s in problem (P) and subsequent problems. In this subsection, we discuss how additional constraint(s) on s, such as a budget constraint, affect our main results reported in Sections 4 through 6.

Case with λ 0
Suppose that the security investments s is required to lie in some convex set S ⊂ R N + in problem (P). Then, the relaxed problems (P R1 ) and (P R2 ) are still convex. However, the pair (s,p) (resp. (s ′ , p ′ )) in Theorem 2 (resp. 3) is a feasible point of problem (P) if and only ifs ∈ S (resp. s ′ ∈ S). For this reason, the convex relaxations (P R1 ) and (P R2 ) are exact ifs and s ′ lie in S and the condition (11) in Lemma 3 holds. In addition, in order to ensure that s (t+1) belongs to S, line 6 of Algorithm 2 needs to be modified as follows: where P S [·] denotes the projection operator onto S.

Case with λ = 0
The finding in Theorem 5 continues to hold when the minimum element s min of the set S exists, with the minimum element s min being the unique optimal point of the problem in (28) when ρ diag(α • s min + δ) −1 B ≤ 1. Hence, when the recovery rates δ are sufficiently large, only the minimum investments given by s min are needed. Obviously, when S = R N + , the minimum element is 0. Also, for a general constraint set S, the problem in (33) is not guaranteed to be feasible, i.e., ρ(diag(α • s + δ) −1 B) > 1 for all s ∈ S. This means that p se (s) > 0 for all s ∈ S and our methods in Sections 4 and 5 can be applied directly to the problem in (28).

Relaxation of Irreducibility of B
Although we suspect that irreducibility of matrix B is a reasonable assumption for many systems of interest, such as enterprise intranets, some systems may not satisfy this assumption. For this reason, here we discuss how relaxing this assumption affects our results.
Note that the irreducibility of B is used to (i) ensure the existence of a unique equilibrium p * (s) ∈ (0, 1) N of (1) as shown in Theorem 1, and (ii) make use of Lemma 2 for our M-matrix based convex formulations.
We relax the assumption that B is irreducible and instead assume that, for every system i ∈ A, either λ i > 0 or there is another system j ∈ A \ {i} with λ j > 0 and a directed path to i in G. Then, the main results in Theorem 1 still hold, i.e., there is a unique equilibrium p * (s) ∈ (0, 1) N of (1) which is strictly positive and can be computed via iteration (2). Moreover, our formulations and results based on exponential cones, including the convex relaxation (P R2 ) and Theorem 3, are still valid because they rely only on the positivity of p(s). However, the convexity of problem (P R1 ) is not guaranteed and requires extending Lemma 2. Finally, when the aforementioned assumption does not hold, the problem is more complicated; our results cannot be applied directly, and it is still an open problem.

CONCLUSIONS
We studied the problem of determining suitable security investments for hardening interdependent component systems of large systems against malicious attacks and infections. Our formulation aims to minimize the average aggregate costs of a system operator based on the steadystate analysis. We showed that the resulting optimization problem is nonconvex, and proposed a set of algorithms for finding a good solution; two approaches are based on convex relaxations, and the other two look for a local minimizer based on RGM and SCP. In addition, we derived a sufficient condition under which the convex relaxations are exact. Finally, we evaluated the proposed algorithms and demonstrated that, although the original problem is nonconvex, local minimizers found by the RGM and SCP methods are good solutions with only small optimality gaps. In addition, as predicted by our analytical results, when the infection costs are high, the optimal points of convex relaxations solve the original nonconvex problem.

APPENDIX A PROOF OF THEOREM 1
The existence and uniqueness of an equilibrium can be established following the same approach in [18], and we omit its proof here.
We now show that p * can be obtained using the iteration in (2) of the theorem. Given λ 0, δ, s ≥ 0, it follows that the steady state p * satisfies (5), which can be rewritten as with u = α • s + δ. This is equivalent to This suggests that p * is a fixed point of the map H. It can be verified that H satisfies the following properties of a standard interference function [41]: for any p ≥ p ′ ≥ 0, (i) H(p) > 0, (ii) H(p) ≥ H(p ′ ), and (iii) cH(p) > H(cp) for c > 1. As a result, the unique fixed point of H can be computed using the iteration starting from any p 0 ≥ 0. The convergence of the synchronous version (as well as asynchronous one) can be found in [41]. However, little is known about the convergence rate except for when a stronger condition than (iii) is used. For example, linear convergence is achieved when H is a contraction [9]. Here, we demonstrate linear convergence using a different, more direct argument. First, for any p ≥ p ′ ≥ p * , we have where Ξ = diag( 1−p * u ), and the second equality is a consequence of 1 − p * = u u+λ+Bp * from (41). Here, for a vector x, we use Next, we show that ρ(ΞB) ≤ 1, which then implies To this end, note ρ(ΞB) ≤ ρ Ξ diag( λ p * ) + B = 1, where the inequality follows from the monotonicity of the spectral radius of nonnegative matrices [13,Thm 8.1.18], and the equality from the fact that p * is a strictly positive eigenvector of Ξ diag( λ p * ) + B corresponding to eigenvalue 1, i.e., = p * .

APPENDIX B A PROOF OF CONVEXITY OF Ω
We will prove this convexity claim by using the convexity of the dominant eigenvalue of an essentially nonnegative matrix (also known as a Metzler matrix, i.e., off-diagonal elements are nonnegative). In particular, it is known that any Metzler matrix A ∈ R n×n has an eigenvalueσ(A) called the dominant eigenvalue that is real and greater than or equal to the real part of any other eigenvalue of A. The following result was shown in [6].
Theorem [6]: Let A be a Metzler matrix and D be a diagonal matrix. Thenσ(A + D) is a convex function of D.
This theorem implies that σ(diag(z) − B) is concave in z.

APPENDIX C FORMULATION OF (12) AS A MATRIX BALANCING PROBLEM
First, recall that a matrix A ∈ R n×n + is balanced if A1 = A T 1. A positive diagonal matrix X is said to balance a matrix A ∈ R n×n + if XAX −1 is balanced. The following lemma holds [39]: Lemma 6. For any A ∈ R n×n + , let G A denote the weighted directed graph associated with A. Then, the following statements are equivalent: (i) A positive diagonal matrix X that balances A exists. (ii) The graph G A is strongly connected. (iii) There exists a positive vector x that minimizes the sum 1≤i,j≤n a ij x i x −1 j . Moreover, matrix X and vector x are unique up to scalars.
We now proceed to prove that (12) is a matrix balancing problem. Suppose that σ(diag(z) − B) = 0, which is equivalent to the condition that diag(z) − B is a singular M-matrix (i.e., z ∈ ∂Ω). Then, there existsB ∈ R N ×N + such that diag(z) − B = ρ(B)I −B. SinceB is also irreducible, it follows from the Perron-Frobenius theorem [13] that there is x > 0 such thatBx = ρ(B)x. As a result, Therefore, we can write (12) as follows.
which, by Lemma 6, is the problem of balancing diag(h)B.

APPENDIX D A PROOF OF THEOREM 3
The first inequality in (17) is obvious because (P R2 ) is a relaxation of (P). The second inequality follows from the feasibility of (s ′ , p ′ ) as shown below.
To show that the above equality holds, we will prove that the equalities hold, i.e., the inequality constraints of t and U are active at the solution x + R of problem (P R2 ). Let µ s , µ p , µ p , µ t , µ y ∈ R N + and Φ ∈ R N ×N + be the Lagrangian multipliers associated with inequality constraints, and σ ∈ R N with the equality constraint in (14) of (P R2 ). Then, the Lagrangian function is given by Since the problem is convex, the following Karush-Kuhn-Tucker (KKT) conditions are necessary and sufficient for optimality.

D.2 Equivalence of Relaxed Problems (P R1 ) and (P R2 )
We will show that (P R1 ) and (P R2 ) are equivalent. Recall Let v := (diag(z) − B) −1 λ. Then, we have Since B is irreducible, Lemma 1-(f) implies that if v > 0, we have z ∈ Ω because λ 0. In addition, we can show that if z ∈ Ω, then (diag(z) − B) −1 > 0 and thus v > 0. 3 Thus, z ∈ Ω if and only if v > 0 and, when v > 0, (48) is equivalent to Rearranging the terms, we get Substituting the expressions for v and z in the first and second constraint, respectively, and replacing the constraint z ∈ Ω with v > 0, we can reformulate (P R1 ) as follows.
Next, we will show that this problem is equivalent to (P R2 ). Recall that, as we showed in the previous subsection, at any optimal point x + = (s + , p + , y + , t + , U + ) of (P R2 ), equalities in (18) hold. Thus, (P R2 ) can be reduced to the following problem after eliminating t and U .
Clearly, (51) is identical to (50) after replacing e −y with v. This proves the equivalence of (P R1 ) and (P R2 ). 3. We can show that the inverse is strictly positive as follows: z ∈ Ω implies that diag(z) − B = µ(I − A) for some µ > 0 and A ∈ R N×N + , which is irreducible with positive diagonal elements and satisfies

APPENDIX E A PROOF OF THEOREM 4
In order to prove the theorem, it suffices to show that at any optimal point x + = (s + , y + , t + , U + ) of the relaxed problem in (27), the constraints (27b) and (27c) are active. In other words, t + = λ • e y + , U + = diag(e y + )Bdiag(e −y + ). (52) The proof is similar to that in Appendix D. Let µ s , µ t ∈ R N + and Φ ∈ R N ×N + be the Lagrange multipliers associated with inequality constraints (27b) and (27c), and σ ∈ R N those associated with the equality constraint (27a) of the relaxed problem. The Lagrangian is given by Since the problem is convex, the necessary and sufficient KKT conditions are given by the following: From (54a) and (54b), By combining this and (54c), we obtain which is strictly positive because c i > 0. Thus, σ > 0, which together with (55) and the slackness conditions in (54d) and (54e), implies that (52) must hold.

APPENDIX F PROOF OF THEOREM 6
Suppose (s L , p L , y L , U L ) is an optimal point of (P R3 ), and let C L = w(s L ). There are two possibilities: C L = C or C L < C. We shall consider them separately below.
• C L = C: In this case, we have where the inequality follows from the assumption C < C * and, hence, p se (s L ) > 0. Since the relaxed problem (P R3 ) is convex, f L (C) is a convex function of C [5]. Thus, we have Combining this with (56) yields Suppose that s 0 is an optimal point of (33). Then, (57) implies that s 0 is O(ǫ)-suboptimal for the original problem in (28).
• C L < C: Since the constraint w(s L ) ≤ C is inactive, it follows that f L (C L ) = f L (C ′ ) for all C ′ ∈ [C L , C]. Because C can be made arbitrarily close to C * , this also tells us that f L (C L ) = f * L , providing us with a following lower bound: This lower bound, together with an upper bound we can obtain as explained in subsection 6.2.2, can be used to quantify how close to optimal a feasible solution is.

APPENDIX G A PROOF OF THEOREM 7
We prove the first part by contradiction: suppose that there exists some C < C * such that, at an optimal point (s L , p L , y L , U L ), we have C L < C. Consider for some γ > 0 such that w(s ′ ) < C. We can find such positive γ because w is continuous and increasing, and w(s L ) = C L < C. We can verify that the tuple (s ′ , p ′ , y ′ , U L ) is a feasible point of the relaxed problem (P R3 ). Moreover, The second inequality holds because γ > 0, p L > 0, s ′ ∈ S C and the condition (38) is true. But, this contradicts the assumed optimality of (s L , p L , y L , U L ). Hence, we have C L = C.
To prove the second part of theorem, assume that (38) holds for all s ∈ S C * . Take any sequence of positive ǫ k , k ∈ N, which satisfies lim k→∞ ǫ k = 0, and consider the sequence {C k = C * − ǫ k : k ∈ N}. Because S C k ⊂ S C * , from (37), we obtain Since the right-hand side goes to 0 as k → ∞, f * 0 = C * and s 0 is an optimal point of the original problem in (28).

APPENDIX H REFORMULATION OF PROBLEM (33) AS A CONVEX (EXPONENTIAL CONE) PROBLEM
Recall that the optimization problem is given by By Perron-Frobenius theory [13], the above spectral radius constraint is equivalent to the existence of a positive vector x such that diag(α • s + δ) −1 Bx ≤ x. With a change of variable y = log(x), problem (P C ) can be expressed as follows.
This is a convex problem and can be turned into an exponential cone optimization problem as follows with U = [u ij ] ∈ R N ×N : min s≥0,y,U w(s) (60) s.t. U 1 ≤ α • s + δ (u ij , 1, y i − y j + log b ij ) ∈ K exp , ∀(i, j) ∈ E This problem can be solved efficiently using convex solvers [1], [8], provided that N is not large. Below, we describe a first-order algorithm for solving the problem in (60) based on our previous matrix balancing algorithm. This algorithm can deal with large-scale problems in a parallel and distributed fashion. For our discussion, without loss of generality, we assume α = 1. Let us consider the Lagrangian given by L(y, s, θ) = w(s) + is a convex problem. Suppose that the Slater's condition holds for problem (58) so that we have strong duality. Then, we can solve (58) by solving the dual problem instead, for example, using a first-order method with subgradients.

APPENDIX I A PROOF OF THE CONTINUITY OF STABLE EQUILIB-RIUM p se
We will prove the claim by contradiction. Suppose that p se is not continuous at some point s ′ ∈ R N + . Then, there exists a sequence s k , k ∈ N, such that s k → s ′ as k → ∞, but either (a) lim k→∞ p se (s k ) does not exist or (b) the limit exists but is not equal to p se (s ′ ). When the limit does not exist, we can work with a convergent subsequence instead; the existence of a convergent subsequence of p se (s k ), k ∈ N, is guaranteed by the Bolzano-Weierstrass Theorem [38]. For this reason, without loss of generality, we assume lim k→∞ p se (s k ) exists and denote it by p † .
We consider the following two cases separately. Note that p se (s) cannot have both positive and zero entries.
Since g(s ′ , p se (s ′ )) = 0, the implicit function theorem [12,Theorem 3] tells us that, for any set P ′ containing p se (s ′ ), we can find a set S ′ containing s ′ and a function ψ : S ′ → P ′ such that (a) ψ(s ′ ) = p se (s ′ ), (b) g(s, ψ(s)) = 0 for all s ∈ S ′ , and (c) ψ is continuous at s ′ . Since there is at most one positive stable equilibrium satisfying g(s, p) = 0, this implies ψ(s) = p se (s) > 0 for all s in a sufficiently small neighborhood around s ′ and, hence, p † = p se (s ′ ), which is a contradiction.
Case 2: p se (s ′ ) = 0 -In this case, since there is a unique solution to g(s ′ , p) = 0, p † must be equal to 0 from the continuity of g, which contradicts the earlier assumption that p † = p se (s ′ ) = 0.

APPENDIX J A DESCRIPTION OF THE SETUP FOR NUMERICAL STUDIES
The setup for running our algorithms is as follows.
• M-matrix + OPTI : In line 3 of Algorithm 1 for solving (P R1 ) with Ω replaced byΩ(z (t) ), we use a general interior point optimizer from package [8] with relative convergence tolerance set to 10 −5 and Hessian matrices approximated by a quasi-Newton algorithm. The solver initial pointx  R according to (7). When computingp (0) R , the iteration in (2) is run until either ||p k+1 − p k || / ||p k || ≤ 10 −7 or k = 500. In addition, we approximate the set of active constraints ofz R (line 4 of Algorithm 1) using I ac = {i ∈ A | [z R ] i − [z] i ≤ 10 −3 }, and selecth = 10. It is important to note that the OPTI package [8] does not exploit/support multithreading; this is also one of the disadvantages compared to MOSEK.
• Exp-cone + MOSEK : We set the relative gap tolerance of the the interior point optimizer in MOSEK to be 10 −5 and do not use MOSEK's presolve procedure as it is timeconsuming for large networks. Moreover, by default, the interior-point optimizer in MOSEK is parallelized and automatically exploits a maximum number of threads.
4. If (s,p) = (s(x * R ),p(x * R )) given in Theorem 2 is available, it can be used as an initial point. Here, we choose (0, p * (0)) to be the initial point for numerical comparisons.
We employ the backtracking line search algorithm with γ 0 = 0.5, a shrinking factor of 0.85, and the Armijo condition parameter set to 10 −4 .
• K-Exp SCP : We note that for conic optimization problems, MOSEK currently offers only an interior-point type optimizer that cannot take advantage of a previous optimal solution [1], requiring a cold-start at each outer iteration. The relative gap tolerance of the interior point optimizer in MOSEK is set to be 10 −5 , and the presolve procedure was ignored. The relative tolerance errors for the outer updates are such that s (t+1) −s (t) s (t) ≤ 10 −6 or |F (s (t+1) )−F (s (t) )| F (s (t) ) ≤ 10 −6 .
• M-matrix SCP : We use Algorithm 3 in combination with the OPTI package. We approximate Ω (t) by Ω (t) mb = s ∈ R n + |s ≥ max{s, δ + λ + Bp (t) } , where s = arg min s∈R n which can be converted into a matrix balancing problem. Here, we set the relative gap tolerance of the interior point optimizer in OPTI to be 10 −5 and the outer tolerance to be the same as in K-Exp SCP. We use a warm-start at each outer iteration, except for the first one.
• MATLAB fmincon: We attempted to use the built-in fmincon function in MATLAB to solve (P) directly in both variables (s, p), but found it inefficient compared to the following form: where s(p) = (p −1 − 1) • (λ + Bp)− δ, obtained from (4b).
We used sqp and interior-point algorithms to solve this problem with constraint and optimality tolerances set to 10 −5 .