key: cord-0146082-cmh5q1h8
authors: Ye, Lintao; Par'e, Philip E.; Sundaram, Shreyas
title: Parameter Estimation in Epidemic Spread Networks Using Limited Measurements
date: 2021-05-11
journal: nan
DOI: nan
sha: 40390415bfa16c5d1c50eea2db65b5a5268e48d0
doc_id: 146082
cord_uid: cmh5q1h8

We study the problem of estimating the parameters (i.e., infection rate and recovery rate) governing the spread of epidemics in networks. Such parameters are typically estimated by measuring various characteristics (such as the number of infected and recovered individuals) of the infected populations over time. However, these measurements also incur certain costs, depending on the population being tested and the times at which the tests are administered. We thus formulate the epidemic parameter estimation problem as an optimization problem, where the goal is to either minimize the total cost spent on collecting measurements, or to optimize the parameter estimates while remaining within a measurement budget. We show that these problems are NP-hard to solve in general, and then propose approximation algorithms with performance guarantees. We validate our algorithms using numerical examples.

Models of spreading processes over networks have been widely studied by researchers from different fields (e.g., [14, 23, 6, 3, 25, 20] ). The case of epidemics spreading through networked populations has received a particularly significant amount of attention, especially in light of the ongoing COVID-19 pandemic (e.g., [20, 21] ). A canonical example is the networked SIR model, where each node in the network represents a subpopulation or an individual, and can be in one of three states: susceptible (S), infected (I), or recovered (R) [18] . There are two key parameters that govern such models: the infection rate of a given node, and the recovery rate of that node. In the case of a novel virus, these parameters may not be known a priori, and must be identified or estimated from gathered data, including for instance the number of infected and recovered individuals in the network at certain points of time. For instance, in the COVID-19 pandemic, when collecting the data on the number of infected individuals or the number of recovered individuals in the network, one possibility is to perform virus or antibody tests on the individuals, with each test incurring a cost. Therefore, in the problem of parameter estimation in epidemic spread networks, it is important and of practical interest to take the costs of collecting the data (i.e., measurements) into account, which leads to the problem formulations considered in this paper. The goal is to exactly identify (when possible) or estimate the parameters in the networked SIR model using a limited number of measurements. Specifically, we divide our analysis into two scenarios: 1) when the measurements (e.g., the number of infected individuals) can be collected exactly without error; and 2) when only a stochastic measurement can be obtained.

Under the setting when exact measurements of the infected and recovered proportions of the population at certain nodes in the network can be obtained, we formulate the Parameter Identification Measurement Selection (PIMS) problem as minimizing the cost spent on collecting the measurements, while ensuring that the parameters of the SIR model can be uniquely identified (within a certain time interval in the epidemic dynamics). In settings where the measurements are stochastic (thereby precluding exact identification of the parameters), we formulate the Parameter Estimation Measurement Selection (PEMS) problem. The goal is to optimize certain estimation metrics based on the collected measurements, while satisfying the budget on collecting the measurements.

The authors in [22, 34] studied the parameter estimation problem in epidemic spread networks using a "Susceptible-Infected-Susceptible (SIS)" model of epidemics. model. When exact measurements of the infected proportion of the population at each node of the network can be obtained, the authors proposed a sufficient and necessary condition on the set of the collected measurements such that the parameters of the SIS model (i.e., the infection rate and the recovery rate) can be uniquely identified. However, this condition does not pose any constraint on the number of measurements that can be collected.

In [24] , the authors considered a measurement selection problem in the SIR model. Their problem is to perform a limited number of virus tests among the population such that the probability of undetected asymptotic cases is minimized. The transmission of the disease in the SIR model considered in [24] is characterized by a Bernoulli random variable which leads to a Hidden Markov Model for the SIR dynamics.

Finally, our work is also closely related to the sensor placement problem that has been studied for control systems (e.g., [19, 37, 36] ), signal processing (e.g., [7, 35] ), and machine learning (e.g., [17] ). The goal of these problems is to optimize certain (problem-specific) performance metrics of the estimate based on the measurements of the placed sensors, while satisfying the sensor placement budget constraints.

First, we show that the PIMS problem is NP-hard, which precludes polynomial-time algorithms for the PIMS problem (if P = NP). By exploring structural properties of the PIMS problem, we provide a polynomial-time approximation algorithm which returns a solution that is within a certain approximation ratio of the optimal. The approximation ratio depends on the cost structure of the measurements and on the graph structure of the epidemic spread network. Next, we show that the PEMS problem is also NP-hard. In order to provide a polynomial-time approximation algorithm that solves the PEMS problem with performance guarantees, we first show that the PEMS problem can be transformed into the problem of maximizing a set function subject to a knapsack constraint. We then apply a greedy algorithm to the (transformed) PEMS problem, and provide approximation guarantees for the greedy algorithm. Our analysis for the greedy algorithm also generalizes the results from the literature for maximizing a submodular set function under a knapsack constraint to nonsubmodular settings. We use numerical examples to validate the obtained performance bounds of the greedy algorithm, and show that the greedy algorithm performs well in practice.

The sets of integers and real numbers are denoted as Z and R, respectively. For a set S, let |S| denote its cardinality. For any n ∈ Z ≥1 , let [n] {1, 2, . . . , n}. Let 0 m×n denote a zero matrix of dimension m × n; the subscript is dropped if the dimension can be inferred from the context. For a matrix P ∈ R n×n , let P , tr(P ) and det(P ) be its transpose, trace and determinant, respectively. The eigenvalues of P are ordered such that |λ 1 (P )| ≥ · · · ≥ |λ n (P )|. Let P ij (or (P ) ij ) denote the element in the ith row and jth column of P , and let (P ) i denote the ith row of P . A positive semidefinite matrix P ∈ R n×n is denoted by P 0.

Suppose a disease (or virus) is spreading over a directed graph G = {V, E}, where V [n] is the set of n nodes, and E is the set of directed edges (and self loops) that captures the interactions among the nodes in V. Here, each node i ∈ V is considered to be a group (or population) of individuals (e.g., a city or a country). A directed edge from node i to node j, where i = j, is denoted by (i, j). For all i ∈ V, denote N i {j : (j, i) ∈ E} andN i {j : (j, i) ∈ E} ∪ {i}. For all i ∈ V and for all k ∈ Z ≥0 , we let s i [k], x i [k] and r i [k] represent the proportions of population of node i ∈ V that are susceptible, infected and recovered at time k, respectively. To describe the dynamics of the spread of the disease in G, we will use the following discrete-time SIR model (e.g., [11] ):

where β ∈ R ≥0 is the infection rate of the disease, δ ∈ R ≥0 is the recovery rate of the disease, h ∈ R ≥0 is the sampling parameter, and a ij ∈ R ≥0 is the weight associated with edge (j, i). Let A ∈ R n×n be a weight matrix, where A ij = a ij for all i, j ∈ V such that j ∈N i , and

T ∈ R n , for all k ∈ Z ≥0 . Throughout this paper, we assume that the weight matrix A ∈ R n×n and the sampling parameter h ∈ R ≥0 are given. Using similar arguments to those in [11] , one can show that

for all i ∈ V and for all k ∈ Z ≥0 . We also have the following result that characterizes properties of x i [k] and r i [k] in the SIR model over G given by Eq. (1); the proof can be found in Section 7.1 in the Appendix. 

Throughout this section, we assume that S I , S H ⊆ V defined in Definition 3.4 are known, i.e., we know the set of nodes in V that have infected individuals initially.

Given exact measurements of x i [k] and r i [k] for a subset of nodes, our goal is to estimate (or uniquely identify, if possible) the unknown parameters β and δ. Here, we consider the scenario where collecting the measurement of x i [k] (resp., r i [k]) at any node i ∈ V and at any time step k ∈ Z ≥0 incurs a cost, denoted as c k,i ∈ R ≥0 (resp., b k,i ∈ R ≥0 ). Moreover, we can only collect the measurements of x i [k] and r i [k] for k ∈ {t 1 , t 1 + 1, . . . , t 2 }, where t 1 , t 2 ∈ Z ≥0 are given with t 2 > t 1 . Noting that Lemma 3.5 provides a (sufficient and necessary) condition under which x i [k] = 0 (resp., r i [k] = 0) holds, we see that one does not need to collect a measurement of x i [k] (resp., r i [k]) if x i [k] = 0 (resp., r i [k] = 0) from Lemma 3.5. Given time steps t 1 , t 2 ∈ Z ≥0 with t 2 > t 1 , we define a set

which represents the set of all candidate measurements from time step t 1 to time step t 2 . To proceed, we first use Eq. (1b)-(1c) to obtain

. . .

1 Note that for the case when di = 0, i.e., i ∈ SI , part (b) implies xi[k] > 0 for all k ≥ 0.

where Φ x

and Φ r

We can then view Eq. (3) as a set of 2(t 2 − t 1 )n equations in β and δ. Noting that

as argued in Section 3, we see that the coefficients in the set of equations in β and δ given by Eq. (3) 

Next, let I ⊆ I t 1 :t 2 denote a measurement selection strategy, where I t 1 :t 2 is given by Eq. (2). We will then consider identifying β and δ using measurements contained in I ⊆ I t 1 :t 2 . To illustrate our analysis, given any i ∈ V and any k ∈ {t 1 , . . . , t 2 − 1}, we first consider the following equation from Eq. (3):

where

, and we index the equation in Eq. (3) corresponding to Eq. (6) as (k, i, x). Note that in order to use Eq. (6) in identifying β and δ, one needs to determine the coefficients (i.e., the terms other than β and δ) in the equation. Also note that in order to determine the coefficients in equation (k, i, x), one can use the measurements contained in I ⊆ I 

where we index the above equation as (k, i, r In general, let us denote the following two coefficient matrices corresponding to equations (k, i, x) and (k, i, r) in Eq. (3), respectively:

for all k ∈ {t 1 , . . . , t 2 − 1} and for all i ∈ V. Moreover, given any measurement selection strategy I ⊆ I t 1 :t 2 , we let

be the set that contains indices of the equations from Eq. (3) that can be potentially used in identifying β and δ, based on the measurements contained in I. In other words, the coefficients on the left-hand side of equation (k, i, x) (resp., (k, i, r)) can be determined using the measurements from I and using Lemma 3.5, for all (k, i, x) ∈Ī (resp., (k, i, r) ∈Ī). Let us now consider the coefficient matrix Φ x k,i (resp., Φ r k,i ) corresponding to (k, i, x) ∈Ī (resp., (k, i, r) ∈Ī). One can then show that it is possible that there exist equations inĪ whose coefficients cannot be (directly) determined using the measurements contained in I or using Lemma 3.5, where the undetermined coefficients come from the first element in Φ x k,i given by Eq. (8a). Nevertheless, it is also possible that one can perform algebraic operations among the equations inĪ such that the undetermined coefficients get cancelled. Formally, we define the following.

Definition 4.1. Consider a measurement selection strategy I ⊆ I t 1 :t 2 , where I t 1 :t 2 is given by Eq. (2) .

k,i and Φ r k,i are given by (8) andĪ is given by Eq. (9) . The resulting matrix is denoted as Φ(I) ∈ R |Ī|×2 . Moreover, defineΦ(I) to be the set that contains all the matrices Φ ∈ R 2×2 such that (Φ) 1 and (Φ) 2 can be obtained via algebraic operations among the rows in Φ(I), and the elements in (Φ) 1 and (Φ) 2 can be fully determined using the measurements from I ⊆ I t 1 :t 2 and using Lemma 3.5.

In other words, Φ ∈Φ(I) corresponds to two equations (in β and δ) obtained from Eq. (3) such that the coefficients on the right-hand side of the two equations can be determined using the measurements contained in I and using Lemma 3.5 (if the coefficients contain x i [k] = 0 or r i [k] = 0). Moreover, one can show that the coefficients on the left-hand side of the two equations obtained from Eq. (3) corresponding to Φ can also be determined using measurements from I and using Lemma 3.5. Putting the above arguments together, we see that given a measurement selection strategy I ⊆ I t 1 :t 2 , β and δ can be uniquely identified if and only if there exists Φ ∈Φ(I) such that rank(Φ) = 2. Equivalently, denoting

where r max (I) 0 ifΦ(I) = ∅, we see that β and δ can be uniquely identified using the measurements from I ⊆ I t 1 :t 2 if and only if r max (I) = 2.

Remark 4.2. Note that if a measurement selection strategy I ⊆ I t 1 :t 2 satisfies r max (I) = 2, it follows from the above arguments that |Ī| ≥ 2, i.e., Φ(I) ∈ R |Ī|×2 has at least two rows, whereĪ is defined in Eq. (9) .

Recall that collecting the measurement of x i [k] (resp., r i [k]) at any node i ∈ V and at any time step k ∈ Z ≥0 incurs cost c k,i ∈ R ≥0 (resp., b k,i ∈ R ≥0 ). Given any measurement selection strategy I ⊆ I t 1 :t 2 , we denote the cost associated with I as c(I)

We then define the Parameter Identification Measurement Selection (PIMS) problem in the perfect measurement setting as follows. (1) with a directed graph G = {V, E}, a weight matrix A ∈ R n×n , a sampling parameter h ∈ R ≥0 , and sets S I , S H ⊆ V defined in Definition 3.4. Moreover, consider time steps t 1 , t 2 ∈ Z ≥0 with t 1 < t 2 , and a cost c k,i ∈ R ≥0 of measuring x i [k] and a cost b k,i ∈ R ≥0 of measuring r i [k] for all i ∈ V and for all k ∈ {t 1 , . . . , t 2 }. The PIMS problem is to find I ⊆ I t 1 :t 2 that solves min

where I t 1 :t 2 is defined in Eq. (2), c(I) is defined in Eq. (11) , and r max (I) is defined in Eq. (10).

We have the following result; the proof is included in Section 7.2 in the Appendix. Theorem 4.4 indicates that there is no polynomial-time algorithm that solves all instances of the PIMS problem optimally (if P = NP). Moreover, we note from the formulation of the PIMS problem given by Problem 4.3 that for a measurement selection strategy I ⊆ I t 1 :t 2 , one needs to check if max Φ∈Φ(I) rank(Φ) = 2 holds, before the corresponding measurements are collected. However, in general, it is not possible to calculate rank(Φ) when no measurements are collected. In order to bypass these issues, we will explore additional properties of the PIMS problem in the following.

We start with the following result. 

To prove part (a), consider any i 1 ∈ S I and any i 2 ∈ V with d i 2 = ∞, where we note x i 1 [0] > 0 and a i 1 i 1 > 0 from the definition of S I . We then see from Lemma 3.5(a)-(b) that s i 1 [k 1 ] > 0 and

We then prove part (b). Considering any i 1 ∈ S and any i 2 ∈ V with d 2 = ∞, we see from the definition of S that N i 1 = ∅ and there exists j ∈ N i 1 such that d j = ∞. Letting j 1 be a node in N i 1 such that d j 1 = min{d j : j ∈ N i 1 } = ∞, we note from Lemma 3.5(a) that x j 1 [k 1 ] > 0 for all k 1 ≥ min{d j : j ∈ N i 1 }. Also note that a i 1 j 1 > 0 from Assumption 3.2. The rest of the proof of part (b) is then identical to that of part (a).

Recalling the way we index the equations in Eq. (3) (see Eqs. (6)- (7) for examples), we define the set that contains all the indices of the equations in Eq. (3):

Following the arguments in Lemma 4.5, we denote

where S I and S are defined in Lemma 4.5, and d i is defined in Definition 3.4. Next, for all (k, i, x) ∈ Q, we define the set of measurements that are needed to determine the coefficients in equation (k, i, x) (when no other equations are used) to be

where I t 1 :t 2 is defined in Eq. (2) . Similarly, for all (k, i, r) ∈ Q, we define

Moreover, let us denote

Similarly to Eq. (11), denote the sum of the costs of the measurements from I((

Based on the above arguments, we propose an algorithm defined in Algorithm 1 for the PIMS problem. Note that Algorithm 1 finds an equation from Q 1 and an equation from Q 2 such that the sum of the costs of the two equations is minimized, where Q 1 and Q 2 are defined in Eq. (15) and Eq. (16), respectively. Proposition 4.6. Consider an instance of the PIMS problem under Assumptions 3.1-3.2. Algorithm 1 returns a solution I((k 1 , i 1 , x), (k 2 , i 2 , r)) to the PIMS problem that satisfies the constraint in (12) , and the following:

where I is an optimal solution to the PIMS problem, Q 1 is defined in Eq. (15) , and c min min{min

Proof. The feasibility of I((k 1 , i 1 , x), (k 2 , i 2 , r)) follows directly from the definition of Algorithm 1 and Lemma 4.5. We now prove (18) . Consider any equations (k, i, x) ∈ Q 1 and (k, i, r) ∈ Q 2 . We have from Eq. (17) the following:

which implies

Next, since I satisfies r max (I ) = 2, we recall from Remark 4.2 |Ī | ≥ 2, wherē

for any w ∈ V, and cannot all be zero, we see that there exists

This further implies |I | ≥ 3. Using similar arguments, one can show that |I | ≥ 3 holds in general, which implies c(I ) ≥ 3c min . Combining the above arguments leads to (18) .

Finally, note that Q 2 and I t 1 :t 2 can be obtained by calling the Breadth-First-Search (BFS) algorithm (e.g., [8] ) |S I | times with O(|S I |(n + |E|)) total time complexity. Also note that the time

In this section, we assume that the initial condition

Nevertheless, our analysis can potentially be extended to cases where the initial condition l is given by a probability distribution.

Here, we consider the scenario where the measurement of

). Note one can express x i [k] in terms of l and θ [β δ] T using (1b). Hence, given l and θ, we can alternatively write p(

for all i ∈ V and for all k ∈ Z ≥1 . Since the initial conditions are assumed to be known, we drop the dependency of p(x i [k]|l, θ) on l, and denote the pmf ofx i [k] as p(x i [k]|θ) for all i ∈ V and for all k ∈ Z ≥1 . Similarly, given l and θ, we denote the pmf ofr i [k] as p(r i [k]|θ) for all i ∈ V and for all k ∈ Z ≥1 . Note that when collecting measurementx i [k] (resp.,r i [k]) under a limited budget, one possibility is to give virus (resp., antibody) tests to a group of randomly and uniformly sampled individuals of the population at node i ∈ V and at time k ∈ Z ≥1 (e.g., [2] ), where a positive testing result indicates that the tested individual is infected (resp., recovered) at time k (e.g., [1] ). Thus, the obtained random measurementsx i [k] andr i [k] and the corresponding pmfs p(x i [k]|θ) and p(r i [k]|θ) depend on the total number of conducted virus tests and antibody tests at node i and at time k, respectively. Consider any node i ∈ V and any time step k ∈ Z ≥1 , where the total population of i is denoted by N i ∈ Z ≥1 and is assumed to be fixed over time. Suppose we are allowed to choose the number of virus (resp., antibody) tests that will be performed on the (randomly sampled) individuals at node i and at time k. Assume that the cost of performing the virus (resp., antibody) tests is proportional to the number of the tests. For any i ∈ V and for any k ∈ {t 1 , . . . , t 2 }, let

be the set of all possible costs that we can spend on collecting the measurementx i [k], where c k,i ∈ R ≥0 and ζ i ∈ Z ≥1 . Similarly, for any i ∈ V and any k ∈ {t 1 , . . . , t 2 }, let

denote the set of all possible costs that we can spend on collecting the measurementr i [k], where b k,i ∈ R ≥0 and η i ∈ Z ≥1 . For instance, ζc k,i can be viewed as the cost of performing virus tests on (20)). Note that ϕ k,i (resp., ω k,i ) is the cost that we spend on collecting measurementx i [k] (resp.,r i [k]), and ϕ k,i = 0 (resp.,

In contrast with the exact measurement case studied in Section 4, it is not possible to uniquely identify β and δ using measurementsx i [k] andr i [k] which are now random variables. Thus, we will consider estimators of β and δ based on the measurements indicated by a measurement selection strategy. Similarly to Section 4, given time steps t 1 , t 2 ∈ Z ≥1 with t 2 ≥ t 1 , define the set of all candidate measurements as

Recalling C k,i and B k,i defined in Eq. (19) and Eq. (20), respectively, we let µ ∈ Z :

where ζ i , η i ∈ Z ≥1 for all i ∈ V. In other words, a measurement selection µ is defined over the integer lattice Z

so that µ is a vector of dimension |U t 1 :t 2 |, where each element of µ corresponds to an element in U t 1 :t 2 , and is denoted as

for all i ∈ V and for all k ∈ {t 1 , . . . , t 2 }. Thus, for any ϕ k,i ∈ C k,i and for any ω k,i ∈ B k,i , there exists µ ∈ M

, respectively, where we drop the dependencies of the pmfs on c k,i and b k,i for notational simplicity.

To proceed, we consider the scenario where measurements can only be collected under a budget constraint given by B ∈ R ≥0 . Using the above notations, the budget constraint can be expressed as

We then consider estimators of θ = [β δ] T based on any given measurement selection µ ∈ M.

Considering any µ ∈ M, we denote

for all i ∈ V and for all λ ∈ {x, r}. For all i ∈ V and for all λ ∈ {x, r} with

we denote the measurement vector indicated by µ ∈ M as

where

, respectively. We then see from Eq. (25) that y(µ) is a random vector whose pmf is denoted as p(y(µ)|θ, µ). Similarly, the pmf of y(U x i ) (resp.,

Assumption 5.1. For any i ∈ V and for any k 1 ,

are independent of each other. Moreover, for any i, j ∈ V (i = j) and for any

The above assumption ensures that measurements from different nodes or from different time steps are independent, and the measurements of x i [k] and r i [k] are also independent. It then follows from Eq. (25) that the pmf of y(µ) can be written as

where we can further write p(y(U

) for all j ∈ U r . In order to quantify the performance (e.g., precision) of estimators of θ based on µ, we use the Bayesian Cramer-Rao Lower Bound (BCRLB) (e.g., [33] ) associated with µ. In the following, we introduce the BCRLB, and explain why we choose it as a performance metric. First, given any measurement µ ∈ M, let F θ (µ) be the corresponding Fisher information matrix defined as

with the expectation E[·] taken with respect to p(y(µ)|θ, µ). Under Assumption 5.1 and some regularity conditions on the pmfs ofx i [k] andr i [k], Eq. (27) can be written as the following (e.g., [13] ):

Consider any estimatorθ(µ) of θ based on a measurement selection µ ∈ M, and assume that we have a prior pdf of θ = [β δ] T , denoted as p(θ). Under some regularity conditions on the pmfs of

, and p(θ), we have (e.g., [32, 33] ):

where Rθ (µ) ∈ R 2×2 is the error covariance of the estimatorθ(µ), the expectation E[·] is taken with respect to p(y(µ)|θ, µ)p(θ), andC(µ) ∈ R 2×2 is the BCRLB associated with the measurement selection µ. The BCRLB is defined as (e.g., [32, 33] )

where E θ [·] denotes the expectation taken with respect to p(θ), F θ (µ) is given by Eq. (27) , and F p ∈ R 2×2 encodes the prior knowledge of θ as

where the second equality holds under some regularity conditions on p(θ) [32] . Thus, the above arguments motivate us to consider (functions of)C(·) as optimization metrics in the measurement selection problem studied in this section, in order to characterize the estimation performance corresponding to a measurement selection µ ∈ M. In particular, we will consider tr(C(·)) and ln det(C(·)), which are widely used criteria in parameter estimation (e.g., [12] ), and are also known as the Bayesian A-optimality and D-optimality criteria respectively in the context of experimental design (e.g., [27] ). First, considering the optimization metric tr(C(·)), we see from the above arguments that (29) directly implies tr(Rθ (µ) ) ≥ tr(C(µ)) for all estimatorsθ(µ) of θ and for all µ ∈ M. Therefore, a measurement selection µ that minimizes tr(C(µ)) potentially yields a lower value of tr(Rθ (µ) ) for an estimatorθ(µ) of θ. Furthermore, there may exist an estimator θ(µ) that achieves the BCRLB (e.g., [32] ), i.e., tr(C(µ)) provides the minimum value of tr(Rθ (µ) ) that can be possibly achieved by any estimatorθ(µ) of θ, given a measurement selection µ. Similar arguments hold for ln det(C(·)). To proceed, denoting

we define the Parameter Estimation Measurement Selection (PEMS) problem. , for all i ∈ V and for all k ∈ {t 1 , . . . , t 2 }; a budget B ∈ R ≥0 ; and a prior pdf p(θ). Supposex i [k] (resp.,r i [k]) is given by a pmf p(x i [k]|θ, ϕ k,i ) (resp., p(r i [k]|θ, ω k,i )), where ϕ k,i ∈ C k,i (resp., ω k,i ∈ B k,i ). The PEMS problem is to find a measurement selection µ that solves

where M is defined in Eq. Note that F p 0 from (31), and f a (0) = tr(C(0)) = tr((F p ) −1 ) and f d (0) = ln det(C(0)) = ln det((F p ) −1 ) from Eq. (30) . We further assume that F p 0 in the sequel, which implies f (µ) > 0 for all µ ∈ M.

In this section, we consider a measurement model with specific pmfs ofx i [k] andr i [k] (e.g., [4] and [11] ). Nonetheless, our analysis can potentially be extended to other measurement models.

Consider any i ∈ V and any k ∈ {t 1 , . . . , t 2 }. Assume that the total population of node i is fixed over time and is denoted as N i ∈ Z ≥1 . Given any measurement selection µ ∈ M with M defined in Eq. (22), we recall from Section 5.1 that µ(x i [k])c k,i can be viewed as the cost of performing virus tests on µ(x i [k])N x i randomly and uniformly sampled individuals in the population of node i ∈ V,

is the proportion of population at node i and at time k that is infected, and x i [k] ∈ [0, 1) under Assumptions 3.1-3.2 as shown by Lemma 3.5. Thus, a randomly and uniformly sampled individual in the population at node i and at time k will be an infected individual (at time k) with probability x i [k], and will be a non-infected (i.e., susceptible or recovered) individual with probability 1 − x i [k]. Supposing the tests are accurate, 2 we see from the above arguments that the obtained number of individuals that are tested positive, i.e., N ixi [k], is a binomial random variable with parameters

Thus, for any i ∈ V and for any k ∈ {t 1 , . . . , t 2 }, the pmf of

where

Moreover, since the weight matrix A ∈ R n×n and the sampling parameter h ∈ R ≥0 are assumed to be given, we see that given θ = [β δ] T and initial

can be obtained using Eq. (1b) for all i ∈ V and for all k ∈ {t 1 , . . . , t 2 }, where we can view x i [k] as a function in the unknown parameter θ. In other words, given l, θ, µ(x i [k]), N x i and N i , one can obtain the right-hand side of Eq. (34) . Again, we only explicitly express the dependency of the pmf ofx i [k] on θ and µ(x i [k]) in Eq. (34) . Following similar arguments to those above, we assume that for any i ∈ V and for any k ∈ {t 1 , . . . , t 2 }, measurement r i [k] has the following pmf:

where

Similarly, we note that the pmf ofr i [k] given in Eq. (35) reduces to p(r i [k] = 0|θ, µ(r i [k])) = 1 when r i [k] = 0. Considering any measurement selection µ ∈ M and any measurementλ i [k] ∈ U t 1 :t 2 , where λ ∈ {x, r} and U t 1 :t 2 is defined in Eq. (21), we have the following: 

Under the measurement model described above, we show that the PEMS problem is also NP-hard, i.e., there exist instances of the PEMS problem that cannot be solved optimally by any polynomial-time algorithm (if P = NP). The proof of the following result is included in Section 7.3 in the Appendix.

Theorem 5.3. The PEMS problem is NP-hard.

Theorem 5.3 motivates us to consider approximation algorithms for solving the PEMS problem. To begin with, we note that the objective function in the PEMS problem can be viewed as a function defined over an integer lattice. We then have f a : M → R ≥0 and f d : M → R ≥0 , where M is defined in Eq. (22) . First, considering f a : M → R ≥0 , we will define a set function f P a : 2M → R ≥0 , wherē M is a set constructed as

(38) In other words, for any i ∈ V and for any k ∈ {t 1 , . . . , t 2 }, we associate elements (

is then defined as

where for any Y ⊆M, we define µ Y ∈ M such that µ Y ( . Also note that f P a (∅) = 0. Following the arguments leading to (37), we define

is taken with respect to the prior pdf p(θ). Given any θ = [β δ] T , we see from the arguments for (37) that

Similarly, one can obtain E θ 

Next, let the cost of (x i [k], l 1 ) be c k,i , denoted as c(x i [k], l 1 ), for all (x i [k], l 1 ) ∈M, and let the cost of (r i [k], l 2 ) be b k,i , denoted as c(r i [k], l 2 ), for all (r i [k], l 2 ) ∈M, where c k,i ∈ R >0 and b k,i ∈ R >0 are given in the instance of the PEMS problem. Setting the cost structure of the elements inM in this way, we establish an equivalence between the cost of a subset Y ⊆M and the cost of µ Y ∈ M, where µ Y is defined above. Similarly, considering the objective function f d : M → R ≥0 in the PEMS problem, we define a set function f P d : 2M → R ≥0 as

, l 2 ) ∈ Y}| for all i ∈ V and for all k ∈ {t 1 , . . . , t 2 }. Note that given an instance of the PEMS problem in Problem 5.2, we can construct the setM with the associated costs of the elements inM in O(n(t 2 − t 1 + 1)(ζ + η)) time, where n is the number of nodes in graph G = {V, E}, and ζ, η ∈ Z ≥1 with ζ i ≤ ζ and η i ≤ η for all i ∈ V. Assuming that ζ and η are (fixed) constants, the construction of the setM with the associated costs takes O(n(t 2 − t 1 + 1)) time, which is polynomial in the parameters of the PEMS problem (Problem 5.2). Based on the above arguments, we further consider the following problem:

where f P (·) can be either of f P a (·) or f P d (·) with f P a (·) and f P d (·) given by in (41) and (42), respectively, and c(Y) y∈Y c(y) for all Y ⊆M. By the way we construct f P (·) and the costs of elements inM, one can verify that Y a ⊆M (resp., Y d ⊆M) is an optimal solution to Problem (P) with f P (·) = f P a (·) (resp., f P (·) = f P d (·)) if and only if µ Y a (resp., µ Y d ) defined above is an optimal solution to (33) in Problem 5.2 with f (·) = f a (·) (resp., f (·) = f d (·)). Thus, given a PEMS instance we can first constructM with the associated cost for each element inM, and then solve Problem (P).

Note that Problem (P) can be viewed as a problem of maximizing a set function subject to a knapsack constraint, and greedy algorithms have been proposed to solve this problem with performance guarantees when the objective function is monotone nondecreasing and submodular 3 (e.g., [16] and [29] ). Before we formally introduce the greedy algorithm for the PEMS problem, we first note from (40)-(42) that given a prior pdf of θ and any Y ⊆M, one has to take the expectation E θ [·] in order to obtain the value of f P (Y). However, it is in general intractable to explicitly calculate the integration corresponding to E θ [·]. Hence, one may alternatively evaluate the value of f P (Y) using numerical integration with respect to θ = [β δ] T (e.g., [28] ). Specifically, a typical numerical integration method (e.g., the trapezoid rule) approximates the integral of a function (over an interval) based on a summation of (weighted) function values evaluated at certain points within the integration interval, which incurs an approximation error (see e.g., [28] , for more details). We then see from (40)-(42) that in order to apply the numerical integration method described above to f P (Y), one has to obtain the values of

∂θ , and ∂r i [k] ∂θ for a given θ (within the integration interval), where i ∈ V and t 1 ≤ k ≤ t 2 with t 1 , t 2 given in an instance of the PEMS problem. Recall that the initial conditions s[0], x[0] and r[0] are assumed to be known. We first observe that for any given θ, the values of x i [k] and r i [k] for all i ∈ V and for all k ∈ {t 1 , . . . , t 2 } can be obtained using the recursions in (1) in O((t 2 − t 1 + 1)n 2 ) time. Next ∂δ ] T , we take the derivative with respect to β on both sides of the equations in (1) and obtain

Considering any given β, we can then use the recursion in (1) ∂δ and ∂r i [k] ∂δ for all i ∈ V and for all k ∈ {t 1 , . . . , t 2 } in O((t 2 − t 1 + 1)n 2 ) time.

Putting the above arguments together and considering the prior pdf of θ, i.e., p(θ), we see from (40)-(42) that for all Y ⊆M, an approximation of f P (Y), denoted asf P (Y), can be obtained in O(n I (t 2 − t 1 + 1)n 2 ) time, where n I ∈ Z ≥1 is the number of points used for the numerical integration method with respective to θ, as we described above. 4 Furthermore, in the sequel, we may assume 5 We are now ready to introduce the greedy algorithm given in Algorithm 2 to solve the PEMS problem, wherê f P (·) ∈ {f P a (·),f P d (·)} andf P (·) is the approximation of f P (·) as we described above. From the definition of Algorithm 2, we see that the number of function calls off P (·) required in the algorithm is O(|M| 2 ), and thus the overall time complexity of Algorithm 2 is given by O(n I (t 2 − t 1 + 1)n 2 |M| 2 ).

We proceed to analyze the performance of Algorithm 2 when applied to the PEMS problem. First, one can observe that f P d (Y) = ln det(F p + H(Y)) − ln det(F p ) in Problem (P) shares a similar form to the one in [30] . Thus, using similar arguments to those in [30] , one can show that f P d (·) is monotone nondecreasing and submodular with f P d (∅) = 0. Noting the assumption that |f P d (Y) − f P d (Y)| ≤ ε/2 for all Y ⊆M, one can show that y given by line 6 of Algorithm 2

for all y ∈ C. Similarly, one can show that max y∈M f P d (y) ≤ f P d (Y 1 ) + ε, where Y 1 is given by line 3 in Algorithm 2. One can then use similar arguments to those for Theorem 1 in [16] and obtain the following result; the detailed proof is omitted for conciseness.

Theorem 5.4. Consider Problem (P) with the objective function f P d : 2M → R ≥0 given by (42). Then Algorithm 2 yields a solution, denoted as Y g d , to Problem (P) that satisfies

where Y d ⊆M is an optimal solution to Problem (P), c min = min y∈M c(y), 6 and ε ∈ R ≥0 satisfies

Algorithm 2 Greedy algorithm for PEMS 1: Input: An instance of PEMS transformed into the form in (P) 2: Output: Y g 3: Find Y 1 ∈ arg max y∈Mf P (y) 4: Initialize Y 2 = ∅ and C =M 5: while C = ∅ do 6:

Find y ∈ arg max y∈Cf

if c(y ) + c(Y 2 ) ≤ B then 8:

In contrast to f P d (·), the objective function f P a (·) is not submodular in general (e.g., [17] ). In fact, one can construct examples where the objective function f P a (Y) = tr((F p ) −1 − tr (F p + H(Y)) −1 ) in the PEMS problem is not submodular. Hence, in order to provide performance guarantees of the greedy algorithm when applied to Problem (P) with f (·) = f P a (·), we will extend the analysis in [16] to nonsubmodular settings. To proceed, note that for all A ⊆ B ⊆M, we have

. Therefore, the objective function f P a (·) is monotone nondecreasing with f P a (∅) = 0. We then characterize how close f P a (·) is to being submodular by introducing the following definition.

Suppose Algorithm 2 is applied to solve Problem (P). For all j ∈ {1, . . . , |Y 2 |}, let Y j 2 = {y 1 , . . . , y j } denote the set that contains the first j elements added to set Y 2 in Algorithm 2, and let Y 0 2 = ∅. The type-1 greedy submodularity ratio of f P a (·) is defined to be the largest γ 1 ∈ R that satisfies

for all A ⊆M and for all j ∈ {0, . . . , |Y 2 |}. The type-2 greedy submodularity ratio of f P a (·) is defined to be the largest γ 2 ∈ R that satisfies

for all j ∈ {0, . . . , |Y 2 |} and for all y ∈M \ Y j 2 such that c(y) + c(Y j 2 ) > B, where Y 1 ∈ arg max y∈Mf P a (y).

Remark 5.6. Note that f P a (·) is monotone nondecreasing as argued above. Noting the definition of γ 1 in (45), one can use similar arguments to those in [5] and show that γ 1 ∈ [0, 1]; if f P a (·) is submodular, then γ 1 = 1. Similarly, one can show that γ 2 ≥ 0. Supposing that Y 1 ∈ arg max y∈M f P a (y), one can further show that if f P a (·) is submodular, then γ 2 ≥ 1.

Note that since we approximate f P a (·) usingf P a (·) as we argued above, we may not be able to obtain the exact values of γ 1 and γ 2 from Definition 5.5. Moreover, finding γ 1 may require an exponential number of function calls of f P a (·) (orf P a (·)). Nonetheless, it will be clear from our analysis below that obtaining lower bounds on γ 1 and γ 2 suffices. Here, we describe how we obtain a lower bound on γ 2 usingf P a (·), and defer our analysis for lower bounding γ 1 to the end of this section, which requires more careful analysis. Similarly to (46), letγ 2 denote the largest real number that satisfiesf

for all j ∈ {0, . . . , |Y 2 |} and for all y ∈M \ Y j 2 such that c(y) + c(Y j 2 ) > B. Noting our assumption that |f P a (Y) − f P a (Y)| ≤ ε/2 for all Y ⊆M (withf P a (∅) = 0), one can see thatγ 2 given by (47) also satisfies (46), which impliesγ 2 ≤ γ 2 . Given Y j 2 for all j ∈ {0, . . . , |Y 2 |} from Algorithm 2,γ 2 can now be obtained via O(|M| 2 ) function calls off P a (·).

Based on Definition 5.5, the following result extends the analysis in [15, 16] , and characterizes the performance guarantees of Algorithm 2 for Problem (P) with f P (·) = f P a (·).

Theorem 5.7. Consider Problem (P) with the objective function f P a : 2M → R ≥0 given by (39). Then Algorithm 2 yields a solution, denoted as Y g a , to Problem (P) that satisfies

where Y a ⊆M is an optimal solution to Problem (P), γ 1 ∈ R ≥0 and γ 2 ∈ R ≥0 are defined in Definition 5.5, c min = min y∈M c(y), c max = max y∈M c(y), and ε ∈ R ≥0 satisfies |f P a (Y) − f P a (Y)| ≤ ε/2 for all Y ⊆M.

Proof. Noting that (48) holds trivially if γ 1 = 0 or γ 2 = 0, we assume that γ 1 > 0 and γ 2 > 0. In this proof, we drop the subscript of f P a (·) (resp.,f P a (·)) and denote f (·) (resp.,f (·)) for notational simplicity. First, recall that for all j ∈ {1, . . . , |Y 2 |}, we let Y j 2 = {y 1 , . . . , y j } denote the set that contains the first j elements added to set Y 2 in Algorithm 2, and let Y 0 2 = ∅. Now, let j l be the first index in {1, . . . , |Y 2 |} such that a candidate element y ∈ arg max y∈Cf

for Y 2 (given in line 6 of Algorithm 2) cannot be added to Y 2 due to c(y ) + c(Y j l 2 ) > B. In other words, for all j ∈ {0, . . . , j l − 1}, any candidate element y ∈ arg max y∈Cf

for Y 2 satisfies c(y ) + c(Y j 2 ) ≤ B and can be added to Y 2 in Algorithm 2. Noting that |f P (Y) − f P (Y)| ≤ ε/2 for all Y ⊆M, one can then show that the following hold for all j ∈ {0, . . . , j l − 1}:

Now, considering any j ∈ {0, . . . , j l − 1}, we have the following:

where (50) follows from the definition of γ 1 in (45), and (51) follows from (49). To obtain (52), we use the fact c(Y a ) ≤ B. Similarly, we obtain (53). Noting that f (·) is monotone nondecreasing, one can further obtain from (53) that

To proceed, let y ∈ arg max y∈Cf

be the (first) candidate element for Y 2 that cannot be added to Y 2 due to c(y ) + c(Y j l 2 ) > B, as we argued above. Similarly to (49), one can see that

holds for all y ∈M \ Y j l 2 . LettingȲ j l +1 2 {y } ∪ Y j l 2 and following similar arguments to those leading up to (54), we have

Denoting

, we obtain from (54) and (55) the following:

Unrolling (56) yields

To obtain (57), we use the facts that (1 − c(y j )γ 1 B ) ≤ 1 for all j ∈ [j l + 1] and (1 − c(y )γ 1 B ) ≤ 1, since γ 1 ∈ (0, 1] as we argued in Remark 5.6. To obtain (58), we first note from the way we defined j l that j l + 1 ≤ c(Ȳ j l +1 2 )/c min ≤ (B + c max )/c min . Also noting that |Y a | ≤ B/c min , we then obtain (58).

(e.g., [15] ). We then have from (58) the following:

where (59) follows from c(Ȳ j l +1 2 ) > B. To proceed with the proof of the theorem, we note from the definition of γ 2 in Definition 5.5

whereB B + c max . Since f (·) is monotone nondecreasing, we obtain from (60)

We will then split our analysis into two cases. First, supposing that γ 2 ≥ 1, we see from (61) that at least one of f (

holds. Second, supposing γ 2 < 1 and using similar arguments, we have that at least one off (

Combining the above arguments together, we obtain (48).

Also note that γ 2 ≥ 1 can hold when the objective function f P a (·) is not submodular, as we will see later in our numerical examples.

Remark 5.9. The authors in [31] also extended the analysis of Algorithm 2 to nonsubmodular settings, when the objective function can be exactly evaluated (i.e., ε = 0). They obtained a performance guarantee for Algorithm 2 that depends on a submodularity ratio defined in a different manner. One can show that the submodularity ratios defined in Definition 5.5 are lower bounded by the one defined in [31] , which further implies that the performance bound (when ε is 0) for Algorithm 2 given in Theorem 5.7 is tighter than that provided in [31] .

Finally, we aim to provide a lower bound on γ 1 that can be computed in polynomial time. The lower bounds on γ 1 and γ 2 together with Theorem 5.7 will also provide performance guarantees for the greedy algorithm.

Lemma 5.10. ( [10] ) For any positive semidefinite matrices P, Q ∈ R n×n , λ 1 (P ) ≤ λ 1 (P + Q) ≤ λ 1 (P ) + λ 1 (Q), and λ n (P + Q) ≥ λ n (P ) + λ n (Q).

We have the following result; the proof is included in Section 7.4 in the Appendix.

Lemma 5.11. Consider the set function f P a : 2M → R ≥0 defined in (39). The type-1 greedy submodularity ratio of f P a (·) given by Definition 5.5 satisfies

where Y j 2 contains the first j elements added to Y 2 in Algorithm 2 ∀j ∈ {1, . . . , |Y 2 |} with Y 0 2 = ∅, F p is given by (31) ,

Recalling our arguments at the beginning of this section, we may only obtain approximations of the entries in the (2 by 2) matrix F p + H(Y) for Y ⊆M using, e.g., numerical integration, where H y (resp., F p ) is defined in (40) (resp., (31) ). Specifically, for all Y ⊆M, letĤ(Y) = (F p +H(Y))+E(Y) be the approximation of F p + H(Y), where each entry of E(Y) ∈ R 2×2 represents the approximation error of the corresponding entry in F p + H(Y). Since F p and H(Y) are positive semidefinite matrices, E(Y) is a symmetric matrix. Now, using a standard eigenvalue perturbation result, e.g., Corollary 6.3.8 in [10] , one can obtain that

tr(E(Y) E(Y)) denotes the Frobenius norm of a matrix. It then follows that

where ε ∈ R ≥0 and satisfies E(Y) F ≤ ε for all Y ⊆M. Combining the above arguments with (62) in Lemma 5.11, we obtain a lower bound on γ 1 that can be computed using O(|M| 2 ) function calls ofĤ(·).

Note that one can further obtain from (62):

where z j ∈ arg min y∈M\Y j

. Supposing F p is fixed, we see from (63) that the lower bound on γ 1 would potentially increase if λ 2 (H(z j ))/λ 1 (H(z j )) and λ 2 (H(Y j 2 ))/λ 1 (H(Y j 2 )) increase. Recall that F p given by (31) encodes the prior knowledge that we have about θ = [β δ] T . Moreover, recall from (40) that H(y) depends on the prior pdf p(θ) and the dynamics of the SIR model in (1) . Thus, the lower bound given by Lemma 5.11 and thus the corresponding performance bound of Algorithm 2 given in Theorem 5.7 depend on the prior knowledge that we have on θ = [β δ] T and the dynamics of the SIR model. Also note that the performance bounds given in Theorem 5.7 are worst-case performance bounds for Algorithm 2. Thus, in practice the ratio between a solution returned by the algorithm and an optimal solution can be smaller than the ratio predicted by Theorem 5.7, as we will see in our simulations in next section. Moreover, instances with tighter performance bounds potentially imply better performance of the algorithm when applied to those instances. Similar arguments hold for the performance bounds provided in Theorem 5.4.

To validate the theoretical results in Theorems 5.4 and 5.7, and Lemma 5.11, we consider various PEMS instances. 7 The directed network G = {V, E} is given by Fig. 1(a) . According to the existing literature about the estimated infection and recovery rates for the COVID-19 pandemic (e.g., [26] ), we assume that the infection rate β and the recovery rate δ lie in the intervals [3, 7] and [1, 4] , respectively. Let the prior pdf of β (resp., δ) be a (linearly transformed) Beta distribution with parameters α 1 = 6 and α 2 = 3 (resp., α 1 = 3 and α 2 = 4), where β and δ are also assumed to be independent. The prior pdfs of β and δ are then plotted in Fig. 1(b) and Fig. 1(c) , respectively. We set the sampling parameter h = 0.1. We then randomly generate the weight matrix A ∈ R 5×5 such that First, let us consider PEMS instances with a relatively smaller size. In such instances, we set the time steps t 1 = t 2 = 5, i.e., we only consider collecting measurements at time step t = 5. In the sets C 5,i = {ζc 5,i : ζ ∈ ({0} ∪ [ζ i ])} and B 5,i = {ηb 5,i : η ∈ ({0} ∪ [η i ])}, we let c 5,i = b 5,i and ζ i = η i = 2 for all i ∈ V, and draw c 5,i and b 5,i uniformly randomly from {1, 2, 3}. Here, we can choose to perform 0, 100, or 200 virus (or antibody) tests at a node i ∈ V and at k = 5. In Fig. 2(a) , we consider the objective function f P d (·), given by Eq. (42), in the PEMS instances constructed above, and plot the greedy solutions and the optimal solutions (found by brute force) to the PEMS instances under different values of budget B. Note that for all the simulation results in this section, we obtain the averaged results from 50 randomly generated A matrices as described above, for each value of B. As shown in Theorem 5.4, the greedy algorithm yields a 1 2 (1 − e −1 ) ≈ 0.31 approximation for f P d (·) (in the worst case), and the results in Fig. 2(a) show that the greedy algorithm performs near optimally for the PEMS instances generated above. Similarly, in Fig. 2(b) , we plot the greedy solutions and the optimal solutions to the PEMS instances constructed above under different values of B, when the objective function is f P a (·) given in Eq. (39). Again, the results in Fig. 2(b) show that the greedy algorithm performs well for the constructed PEMS instances. Moreover, according to Lemma 5.11, we plot the lower bound on the submodularity ratio γ 1 of f P a (·) in Fig. 2(c) . Here, we note that the submodularity ratio γ 2 of f P a (·) is always greater than one in the PEMS instances constructed above. Hence, Theorem 5.7 yields a 1 2 (1 − e −γ 1 ) worst-case approximation guarantee for the greedy algorithm, where We then investigate the performance of the greedy algorithm for PEMS instances of a larger size. Since the optimal solution to the PEMS instances cannot be efficiently obtained when the sizes of the instances become large, we obtain the lower bound on the submodularity ratio γ 1 of f P a (·) provided in Lemma 5.11, which can be computed in polynomial time. Different from the smaller instances constructed above, we set t 1 = 1 and t 2 = 5. We let ζ i = η i = 10 for all

where we also set c k,i = b k,i and draw c k,i and b k,i uniformly randomly from {1, 2, 3}, for all k ∈ [5] and for all i ∈ V. Moreover, we modify the parameter of the Beta distribution corresponding to the pdf of β to be α 1 = 8 and α 2 = 3. Here, we can choose to perform 0, 100, 200, ..., or 1000 virus (or antibody) tests at a node i ∈ V and at k ∈ [5] . In Fig. 3(a) , we plot the lower bound on γ 1 obtained from the PEMS instances constructed above. We note that the submodularity ratio γ 2 of f P a (·) is also always greater than one. Hence, Theorem 5.7 yields a 1 2 (1 − e −γ 1 ) worst-case approximation guarantee for the greedy algorithm. We plot in Fig. 3 (b) the approximation guarantee using the lower bound that we obtained on γ 1 .

(a) Bound on γ1 (b) Approx. guarantee of greedy 

We first considered the PIMS problem under the exact measurement setting, and showed that the problem is NP-hard. We then proposed an approximation algorithm that returns a solution to the PIMS problem that is within a certain factor of the optimal one. Next, we studied the PEMS problem under the noisy measurement setting. Again, we showed that the problem is NP-hard. We applied a greedy algorithm to solve the PEMS problem, and provided performance guarantees on the greedy algorithm. We presented numerical examples to validate the obtained performance bounds of the greedy algorithm, and showed that the greedy algorithm performs well in practice.

Fact 3. Consider any i ∈ V and any k 1 ∈ Z ≥0 . If

Let us first prove Fact 1. Consider any i ∈ V and any k ∈ Z ≥0 . Supposing x i [k] > 0, we have from Eq. (1) x

where the first term on the right-hand side of the above equation is positive, since 1 − hδ > 0 from Assumption 3.2, and the second term on the right-hand side of the above equation is nonnegative. It then follows that x i [k + 1] > 0. Repeating the above argument proves Fact 1. We next prove Fact 2. Considering any i ∈ V and any k ∈ Z ≥0 such that x i [k] = 0, we note from Eq. (1) that 

We show that PIMS is NP-hard via a polynomial reduction from the exact cover by 3-sets (X3C) problem which is known to be NP-complete [9] . Problem 7.1. Consider X = {1, 2, . . . , 3m} and a collection Z = {z 1 , z 2 , . . . , z τ } of 3-element subsets of X , where τ ≥ m. The X3C problem is to determine if there is an exact cover for X , i.e., a subcollection Z ⊆ Z such that every element of X occurs in exactly one member of Z . Consider an instance of the X3C problem given by a set X = {1, . . . , 3m} and a collection Z = {z 1 , . . . , z τ } of 3-element subsets of X , where τ ≥ m. We then construct an instance of the PIMS problem as follows. The node set of the graph G = {V, E} is set to be V = {i 0 , i 1 , . . . , i τ } ∪ {j 1 , j 2 , . . . , j 3m }. The edge set of G = {V, E} is set to satisfy that (j q , i l ) ∈ E if q ∈ X is contained in z l ∈ Z, (j q , i 0 ) ∈ E for all q ∈ X , and (i 0 , i 0 ) ∈ E. Note that based on the construction of G = {V, E}, each node i ∈ {i 1 , . . . , i τ } represents a subset from Z in the X3C instance, and each node j ∈ {j 1 , . . . , j 3m } represents an element from X in the X3C instance, where the edges between {i 1 , . . . , i τ } and {j 1 , . . . , j 3m } indicate how the elements in X are included in the subsets in Z. A plot of G = {V, E} is given in Fig. 4 . Accordingly, the weight matrix A ∈ R (3m+τ +1)×(3m+τ +1) is set to satisfy that a i l jq = 1 if q ∈ X is contained in z l ∈ Z, a i 0 jq = 1 for all q ∈ X , and a i 0 i 0 = 1. We set the sampling parameter to be h = 1/(3m + 1). The set S I ⊆ V is set to be S I = V, i.e., x i [0] > 0 for all i ∈ V. We set time steps t 1 = 2 and t 2 = 3. Finally, we set b 2,i = b 3,i = 0 for all i ∈ V, c 2,i l = 1 and c 3,i l = 0 for all l ∈ {1, . . . , τ }, c 2,jq = c 3,jq = m + 1 for all q ∈ X , and c 2,i 0 = c 3,i 0 = 0. Since x i [0] > 0 for all i ∈ V, we see from Lemma 3.5 that x i [k] > 0 and r i [k] > 0 for all i ∈ V and for all k ∈ {2, 3}. Therefore, Lemma 3.5 is no longer useful in determining the coefficients in the equations from Eq. (3).

We claim that an optimal solution, denoted as I , to the constructed PIMS instance satisfies c(I ) ≤ m if and only if the solution to the X3C instance is "yes".

First, suppose the solution to the X3C instance is "yes". Denote an exact cover as Z = {z q 1 , . . . , z qm } ⊆ Z, where {q 1 , . . . , q m } ⊆ {1, . . . , τ }. Let us consider a measurement selection strategy I 0 ⊆ I t 1 :t 2 given by

We then have from Eq.

> 0 for all i ∈ V and for all k ∈ Z ≥0 from Lemma 3.5(a), we consider the following (m + 1) equations from Eq. (3) whose indices are contained inĪ 0 :

where we note N i 0 = {j 1 , . . . , j 3m } from the way we constructed G = {V, E}. Since Z = {z q 1 , . . . , z qm } is an exact cover for X , we see from the construction of G = {V, E} that l∈{1,...,m} N iq l is a union of mutually disjoint (3-element) sets such that l∈{1,...,m} N iq l = {j 1 , . . . , j 3m }. Thus, subtracting the equations in (68) from Eq. (67), we obtain

where we note x i 0 [2] > 0 as argued above. Following Definition 4.1, we stack coefficient matrices

for all l ∈ {1, . . . , m} into a matrix Φ(I 0 ) ∈ R (m+2)×2 , where Φ r k,i and Φ x k,i are defined in (8) . Now, considering the matrix:

s i 0 [2] + l∈{1,...,m}

we see from the above arguments that (Φ 0 ) 1 and (Φ 0 ) 2 can be be obtained via algebraic operations among the rows in Φ(I 0 ), and the elements in (Φ 0 ) 1 and (Φ 0 ) 2 can be determined using the measurements from I 0 . Therefore, we have Φ 0 ∈Φ(I 0 ), whereΦ(I 0 ) is defined in Definition 4.1.

Noting that x i 0 [2] > 0, we have rank(Φ 0 ) = 2, which implies r max (I 0 ) = 2, where r max (I 0 ) is given by Eq. (10) . Thus, I 0 ⊆ I t 1 :t 2 satisfies the constraint in (12) . Since c(I 0 ) = m from the way we set the costs of collecting measurements in the PIMS instance, we have c(I ) ≤ m.

Conversely, suppose the solution to the X3C instance is "no", i.e., for any subcollection Z ⊆ Z that contains m subsets, there exists at least one element in X that is not contained in any subset in Z . We will show that for any measurement selection strategy I ⊆ I t 1 :t 2 that satisfies r max (I) = 2, c(I) > m holds. Equivalently, we will show that for any I ⊆ I t 1 :t 2 with c(I) ≤ m, r max (I) = 2 does not hold. Consider any I ⊆ I t 1 :t 2 such that c(I) ≤ m. Noting that c 2,jq = c 3,jq = m + 1 for all q ∈ X in the constructed PIMS instance, we have x jq [2] / ∈ I and x jq [3] / ∈ I for all q ∈ X . Moreover, we see that I contains at most m measurements from {x i 1 [2] , . . . , x iτ [2]}. To proceed, let us consider any I 1 ⊆ I t 1 :t 2 such that

where {v 1 , . . . , v m } ⊆ {1, . . . , τ }. In other words, I 1 ⊆ I t 1 :t 2 contains m measurements from {x i 1 [2] , . . . , x iτ [2] } and all the other measurements from I t 1 :t 2 that have zero costs. It follows that c(I 1 ) = m. Similarly to (67) and (68), we have the following (m + 1) equations from Eq. (3) whose indices are contained inĪ 1 (given by Eq. (9)):

1 s iv l [2] (x iv l [3] − x iv l [2]) = h w∈N iv l

x w [2] − x iv l [2] s iv l [2] β δ , ∀l ∈ {1, . . . , m}.

Noting that for any subcollection Z ⊆ Z that contains m subsets, there exists at least one element in X that is not contained in any subset in Z as we argued above, we see that there exists at least one element in X that is not contained in any subset in {z v 1 , . . . , z vm }. It then follows from the way we constructed G = {V, E} that there exists w ∈ N i 0 such that w / ∈ N iv l for all l ∈ {1, . . . , m}. Thus, by subtracting the equations in (73) (multiplied by any constants resp.) from Eq. (72), x w [2] will remain on the right-hand side of the equation in (72). Similarly, consider any equation from (73) indexed by (2, i v l , x) ∈Ī 1 , where l ∈ {1, . . . , m}. First, suppose we subtract Eq. (72) multiplied by some positive constant and any equations in (73) other than equation (2, i v l , x) (multiplied by any constants resp.) from equation (2, i v l , x). Since there exists w ∈ N i 0 such that w / ∈ N iv l for all l ∈ {1, . . . , m} as argued above, we see that x w [2] will appear on the right-hand side of equation (2, i v l , x). Next, suppose we subtract any equations in (73) other than equation (2, i v l , x) (multiplied by any constants resp.) from equation (2, i v l , x). One can check that either of the following two cases holds for the resulting equation (2, i v l , x): (a) the coefficients on the right-hand side of equation (2, i v l , x) contain x jq [2] / ∈ I 1 , where q ∈ X ; or (b) the coefficient matrix on the right-hand side of equation (2, i v l , x) is of the form 0 . Again, we stack Φ r k,i ∈ R 1×2 for all (k, i, r) ∈Ī 1 and Φ x k,i ∈ R 1×2 for all (k, i, x) ∈Ī 1 into a matrix Φ(I 1 ), where we note that Φ r k,i is of the form 0 for all (k, i, r) ∈Ī 1 . One can then see from the above arguments that for all Φ ∈ R 2×2 (if they exist) such that (Φ) 1 and (Φ) 2 can be obtained from algebraic operations among the rows in Φ(I 1 ), and the elements in (Φ) 1 and (Φ) 2 can be determined using the measurements from I 1 , rank(Φ) ≤ 1 holds. It follows that r max (I 1 ) < 2, i.e., constraint r max (I 1 ) = 2 in (12) does not hold. Using similar arguments to those above, one can further show that r max (I) < 2 holds for all c(I) ≤ m, completing the proof of the converse direction of the above claim.

Hence, it follows directly from the above arguments that an algorithm for the PIMS problem can also be used to solve the X3C problem. Since X3C is NP-complete, we conclude that the PIMS problem is NP-hard.

We prove the NP-hardness of the PEMS problem via a polynomial reduction from the knapsack problem which is known to be NP-hard [9] . An instance of the knapsack problem is given by a set D = {d 1 , . . . , d τ }, a size s(d) ∈ Z >0 and a value v(d) ∈ Z >0 for each d ∈ D, and K ∈ Z >0 . The knapsack problem is to find D ⊆ D such that d∈D v(d) is maximized while satisfying d∈D s(d) ≤ K. Given any knapsack instance, we construct an instance of the PEMS problem as follows. Let G = {V, E} be a graph that contains a set of n isolated nodes with n = τ and V = [n]. Set the weight matrix to be A = 0 n×n , and set the sampling parameter as h = 1. The time steps t 1 and t 2 are set to be t 1 = t 2 = 1, i.e., only the measurements of x i [1] and r i [1] for all i ∈ V will be considered. 34) and (35), respectively, with N x i = N r i = v(d i ) and N i = max i∈V v(d i ) for all i ∈ V, where Assumption 5.1 is assumed to hold. Finally, let the prior pdf of β ∈ (0, 1) be a Beta distribution with parameters α 1 = 3 and α 2 = 3, and let the prior pdf of δ ∈ (0, 1) also be a Beta distribution with parameters α 1 = 3 and α 2 = 3, where we take β and δ to be independent. Noting that C 1,i = {0, B + 1} in the PEMS instance constructed above, i.e.,x i [k] incurs a cost of B + 1 > B, we only need to consider measurementsr i [1] for all i ∈ V. Moreover, since B 1,i = {0, s(d i )}, a corresponding measurement selection is then given by µ ∈ {0, 1} V . In other words, µ(i) = 1 if measurementr i [1] is collected (with cost s(d i )), and µ(i) = 0 if measurementr i [1] is not collected. We will see that there is a one to one correspondence between a measurementr i [1] in the PEMS instance and an element d i ∈ D in the knapsack instance.

Consider a measurement selection µ ∈ {0, 1} V . Since r i [0] = 0 and x i [0] = 0.5 for all i ∈ V, Eq. (1c) implies r i [1] = 0.5hδ for all i ∈ V, where h = 1. One then obtain from Eq. (28) and Eq. (37) the following: 

Next, noting that β and δ are independent, one can show via Eq. (31) that F p ∈ R 2×2 is diagonal, where one can further show that (F p ) 11 = (F p ) 22 > 0 using the fact that the pdfs of β and δ are Beta distributions with parameters α 1 = 3 and α 2 = 3. Similarly, one can obtain E θ [1/0.5δ(1 − 0.5δ)] > 0. It now follows from Eq. (74) that

where z 1 , z 2 ∈ R >0 are some constants (independent of µ). Note that the objective in the PEMS instance is given by min µ∈{0,1} V f (µ), where f (·) ∈ {f a (·), f d (·)}. First, considering the objective function f a (µ) = tr(C(µ)), whereC(µ) = (E θ [F θ (µ)] + F p ) −1 , we see from Eq. (75) that tr(C(µ)) is minimized (over µ ∈ {0, 1} V ) if and only if i∈supp(µ) N r i µ(i) is maximized. Similar arguments hold for the objective function f d (µ) = ln det(C(µ)). It then follows directly from the above arguments that a measurement selection µ ∈ {0, 1} V is an optimal solution to the PEMS instance if and only if D {d i : i ∈ supp(µ )} is an optimal solution to the knapsack instance. Since the knapsack problem is NP-hard, the PEMS problem is NP-hard.

Protect Purdue-Purdue University's response to COVID-19

Global dynamics of epidemic spread over complex networks

COVID-19 antibody seroprevalence

Guarantees for greedy maximization of non-submodular functions with applications

Epidemic thresholds in real networks

Sparsity-promoting sensor selection for non-linear measurement models

Introduction to algorithms

Computers and intractability

Matrix analysis

A closed-loop framework for inference, prediction and control of SIR epidemics on networks

Sensor selection via convex optimization

Fundamentals of statistical signal processing: Estimation theory

Maximizing the spread of influence through a social network

The budgeted maximum coverage problem

A note on the budgeted maximization of submodular functions

Near-optimal sensor placements in Gaussian processes: Theory, efficient algorithms and empirical studies

On the dynamics of deterministic epidemic propagation over networks

Sensor selection strategies for state estimation in energy constrained wireless sensor networks

Spread of epidemic disease on networks

Analysis and control of epidemics: A survey of spreading processes on complex networks

Analysis, estimation, and validation of discrete-time epidemic processes

Epidemic processes in complex networks

Smart testing and selective quarantine for the control of epidemics

Optimal resource allocation for network protection against spreading processes

The effect of control strategies to reduce social mixing on outcomes of the COVID-19 epidemic in Wuhan, China: a modelling study

Optimal design of experiments

Introduction to numerical analysis

An online algorithm for maximizing submodular functions

On submodularity and controllability in complex dynamical networks

Lqg control and sensing co-design

Van Trees. Detection, estimation, and modulation theory, part I

Van Trees. Detection, estimation, and modulation theory, part IV: Optimum array processing

Overcoming challenges for estimating virus spread dynamics from data

Sensor selection for hypothesis testing: Complexity and greedy algorithms

Resilient sensor placement for Kalman filtering in networked systems: Complexity and algorithms

On the complexity and approximability of optimal sensor selection and attack for Kalman filtering

We first prove part (a). Considering any i ∈ V and any k ∈ Z ≥0 , we note from Eq. (1a) thatfor all i ∈ V as argued in Section 3, and