{'fraction_non_alphanumeric': 0.04475943586938426, 'fraction_numerical': 0.00901372273883987, 'mean_word_length': 4.16750645200934, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Distributed and Recursive Parameter Estimation in Parametrized Linear State-Space Models

We consider a network of sensors deployed to sense a spatio-temporal field and estimate a parameter of interest. We are interested in the case where the temporal process sensed by each sensor can be modeled as a state-space process that is perturbed by random noise and parametrized by an unknown parameter. To estimate the unknown parameter from the measurements that the sensors sequentially collect, we propose a distributed and recursive estimation algorithm, which we refer to as the incremental recursive prediction error algorithm. This algorithm has the distributed property of incremental gradient algorithms and the on-line property of recursive prediction error algorithms. We study the convergence behavior of the algorithm and provide sufficient conditions for its convergence. Our convergence result is rather general and contains as special cases the known convergence results for the incremental versions of the least-mean square algorithm. Finally, we use the algorithm developed in this paper to identify the source of a gas-leak (diffusing source) in a closed warehouse and also report numerical simulations to verify convergence.


I. INTRODUCTION
A sensor network consists of sensors that are spatially deployed to make observations about a process or field of interest. If the process has a temporal variation, the sensors also obtain observations sequentially in time. An important problem in such networks is to use the spatially and temporally diverse measurements In this paper, we consider a parameter estimation problem when each individual sensor observation process can be modeled as a linear state-space process that is parametrized by an unknown parameter of interest, and also perturbed by process and observation noise. State-space models arise directly, or as linear approximations to non-linear models, in many applications. As an example, we will later discuss the problem of localizing a gas-leak in a warehouse.
We propose a distributed and recursive estimation procedure, which is suitable for in-network processing. Each sensor locally processes its own data and shares only a summary in each time slot. The sensors form a cycle and update incrementally, whereby each sensor updates the estimate using its local information and the received estimate from its upstream neighbor, and passes the updated estimate to its downstream neighbor. In this way, there is a reduction in communication within the network at the cost of increased local sensor processing. This can significantly reduce the total network energy used, especially, when the sensors communicate over a wireless medium. Furthermore, the sensor updates are generated recursively from every new measurement using only a summary statistic of the past measurements. This has two benefits. Firstly, the network has a (possibly coarse) estimate at all times, which is important in applications that require the network to react immediately to a stimulus and make decisions on-line. For example, in a network that is deployed to monitor gas leaks the network should raise an alert depending on the level of the leak intensity. An additional benefit is that each sensor can purge its old measurements periodically since only a summary of constant size is used to update the estimates. This can significantly reduce the memory requirements at the sensors.
Our approach is in contrast with traditional estimation methods such as the maximum likelihood and least-squares which are centralized, i.e., the measurements collected by the spatially distributed sensors are routed through the network to a single location (fusion center) where estimates are computed. In this case, the network energy is mainly consumed in routing the measurements to the fusion center, which can be inefficient in terms of energy consumption. The problem of centralized recursive estimation in linear state-space models is an old problem in system identification. We refer the interested reader to [1] for a survey of these methods for linear state-space models. The problem has also generated considerable interest in the neural networks community where the EM algorithm is used as a tool to learn the parameters [2]. A related algorithm is the parallel recursive prediction error algorithm proposed in [3] that updates the components of the parameter vector in parallel.
The literature on distributed estimation is somewhat limited. A distributed maximum-likelihood algorithm is discussed in [4] and a distributed expectation-maximization algorithm is discussed in [5]. In [6], the incremental (sub)gradient algorithms of [7] are used to obtain distributed least-square estimators.
Distributed linear least-squares are discussed in [8] without an explicit point-to-point message routing.
All of these algorithms are distributed but not recursive. This paper extends our earlier work [9], where we considered the problem of recursive and distributed estimation for stationary models. To the best of our knowledge, there is only one other related study [10] that deals with both distributed and recursive estimation. There, incremental versions of the least-mean square algorithm and the recursive least-squares are developed to solve the linear least-squares problem.
In both studies [9] and [10], the models are not auto-regressive.
Our contribution in this paper is the development and convergence analysis of a general distributed recursive algorithm for parameter estimation in parametrized state-space models. Our results are more general than those of [10], which follow as a special case.
The rest of the paper is organized as follows. We formulate the problem in Section II, and then introduce our notation in Section III. We give an overview of the algorithm in Section IV. We then discuss the standard recursive prediction error algorithm algorithm [1] and the incremental gradient algorithm of [7] in Section V. These algorithms are at the heart of our distributed algorithm presented in Section VI, where we also state our main convergence result. We prove the convergence of the algorithm in Appendix A.
We discuss some simple extensions in Section VII. We report some experimental results obtained by our method as employed to localize the sourse in a gas leak problem in Section VIII. We conclude in Section IX.

II. PROBLEM FORMULATION
We consider a network of m sensors, indexed 1, . . . , m, deployed to sense a spatio-temporal diverse field to determine the value of some quantity of interest, denoted by x, with x ∈ ℜ d . We denote the true value of the parameter by x * .
We assume that time is slotted and each sensor sequentially senses the field once in every time slot.
We model the measurement sequence of sensor i as a random process {R i (k; x)} with the following Here, {W i (k; x)} is the process noise, {V i (k)} is the measurement noise, H i is the observation matrix and V i (k + 1) is the measurement noise of sensor i. The process {Θ i (k; x)} can be interpreted as the temporal process obtained by sampling a spatio-temporal diverse field at the location of sensor i. At this point, we do not assume any knowledge on the joint statistics of Θ i (k + 1; x) and Θ j (k + 1; x). 1 We denote by r i (k) the actual measurement collected by sensor i at time slot k, i.e., r i (k) is a The quantities D i (x), H i , Cov(W i (k; x)) and Cov(V i (k)) are available only at sensor i. Moreover, at all sensors a set X is available that satisfies the following properties 1) The set X is closed and convex; 2) The true parameter x * is contained in the set X; 3) The system in (1) is stable, observable and controllable for all x ∈ X.
Note that X may even be the entire ℜ d . The problem is to estimate the parameter x from the collection of sensor measurements {r i (k)} with an algorithm that is: 1) Distributed: Sensor i does not share its raw measurements {r i (k)} with any other sensor.
2) Recursive: At all times, sensor i stores only a summary statistic of a constant size, i.e., the size of the statistic does not increase with the number of measurements collected by the sensor.

III. NOTATION
All the random variables are defined on the same probability space T = (Ω, F, P). If ω ∈ Ω, is the outcome of an experiment, then for a random process {Y (k; x)} that is parametrized by x, we define y(k) =: Y ω (k; x * ), i.e., y(k) is the value of the random variable Y (k; x * ) corresponding to the outcome ω.
According to this notation, r i (k) and θ i (k) are the realizations of R i (k; x * ) and Θ i (k; x * ) that correspond to the same outcome ω.
We let I denote the set of sensors, i.e., I := {1, . . . , m}. Further, we assume that Θ i (k; x) and R i (k; x) are vectors of dimensions q and p, respectively. They are the same for all sensors i ∈ I. 2 We write x)}, which should be viewed as a collection of random variables parametrized by x and not as a function of x. Furthermore, in line with our notation, r k i denotes the realization of R k i (x * ), i.e., r k i denotes the collection {r i (1), . . . , r i (k)}. 1 Even if some information is available we ignore it. This aspect is discussed in detail in Section VII-B. 2 We make this assumption only for the sake of clarity. Our analysis applies to the general case where the dimensions q and p can be sensor dependent.

IV. ALGORITHM OVERVIEW
A standard estimation procedure defines the estimate as the minimum of a suitably defined cost that is a function of the observations and the unknown parameter. For example, the maximum likelihood estimator minimizes the negative of the log-likelihood function. The form of the cost function determines whether there is a distributed and recursive minimization procedure. Further, the cost function also determines other properties of the estimator such as unbiasedness, consistency, minimum variance etc. (see [11]).
Except in some very special estimation problems, it is impossible to find a cost function that supports even a centralized recursive procedure and also generates a 'good' estimate. In this paper, we develop a distributed and recursive estimator that is only consistent, i.e., the estimate converges to the correct value x * as the number of available measurements becomes infinite 3 . The estimates are biased but this is the price that is to be paid to obtain a distributed and recursive procedure. Thus, there are two aspects to the problem. The first is to choose a suitable cost function, and the second is to develop a distributed and recursive minimization procedure. We will first discuss the cost function that is used, and then give an overview of the minimization algorithm.

A. Cost function
Suppose that each sensor has made N measurements and we want to estimate the parameter x from these measurements. As mentioned, the cost is a function of both the available measurements and the unknown parameter. Therefore, we denote the cost function as f N (x; r N ).
For x ∈ X, we assumed that the system in (1) is stable, observable and controllable. The Kalman gain for the system therefore converges to a finite time-invariant value [12]. Let G i (x) be the Kalman gain for the state-space system in (1), which is determined from as the solution to the Riccati equation [1].
with φ i,1 (x; r 0 i ) = µ i (x). Observe that g i,k+1 (x; r k i ) is linear in r k i for each x. Furthermore, for any x ∈ ℜ d , g i,k+1 (x; r k i ) viewed as a function of r k i is an one-step prediction filter (henceforth, referred to as a predictor) for the random process {R i (k + 1; x * )}. Thus, {ĝ i,k+1 (x; r k i )} x∈ℜ d is a predictor family parametrized by x.
We will choose our cost function to be and our estimator to bex We next provide an intuitive explanation as to why this choice of cost function should generate consistent estimates. First, note that the vector g i,k+1 (x * ; r k i ) is the steady-state Kalman predictor for R i (k + 1; x * ) since {r i (k)} is a sample path of the random process {R i (k; x * )}. Among other properties, the steady state Kalman filter is asymptotically optimal in a mean square sense in the class of linear time-invariant predictors. Thus, x * will minimize Here, the expectations are taken with respect to the random process {R k i (x * )}. Therefore, one can expect thatx N , which is the minimum of f N (x; r N ), might in the limit be equal to x * , since x * is the minimum of f (x), the limit of E f N x; R N (x * ) .

B. Overview of the minimization procedure
Observe that f N (x, r N ) can be written as Suppose we are interested in a non-recursive but distributed solution the problem. Then, the incremental gradient algorithm [7] can be used to minimize f N (x; r N ). In each time slot, the algorithm cycles the estimate through the sensor network. Sensor i receives the estimate z i−1,k+1 from sensor i − 1 at time slot k + 1, and generates a new estimate z i,k+1 using ∇f i (z i−1,k+1 ). The new estimate is then passed to sensor i + 1 for i < m, and to sensor 1 for i = m and, thus, the estimate is cycled through the network for each sensor to update. An illustration is shown in Fig. 1. Now consider the complementary problem, i.e., a centralized but recursive solution. By recursive we mean, the algorithm should be able to obtainx N +1 directly fromx N , the new measurements r 1,N +1 , . . . , r m,N +1 , and some summary statistic of the past observations. This is generally not possible.
Nevertheless, it is possible to use the recursive prediction error algorithm of [1] to obtain recursive approximations to {x k } such that the approximate sequence converges to the same limit as {x k }. Thus, while we do not minimize the chosen cost function f N at each step, the new sequence {x k } is still consistent.
The algorithm that we propose is a combination of the incremental gradient algorithm and the recursive prediction error algorithm. We therefore refer to the algorithm developed in this paper as the incremental recursive prediction error algorithm.

V. PRELIMINARIES
We first evaluate some quantities that will be useful to us in the analysis. To make the paper selfcontained, we then discuss the incremental gradient method [7] and the recursive prediction error algorithm [1] in this section.

A. Some more notation
For later reference, we obtain the form of the gradient of the predictor g i,k+1 x; r k i . Define F i (x) = D i (x) − G i (x)H i and for convenience rewrite (2) as follows Let x (ℓ) denote the ℓ-th component of x, and define Thus the gradient By differentiating in (4), we can immediately see that

B. Incremental gradient descent algorithm
For differentiable optimization problem of the form the standard gradient descent method, with projections, generates iterates according to the following rule: Here, the scalar α k+1 > 0 is the step-size and P X denotes the projection onto the set X. This method is centralized in the sense that it requires the gradient information of each f i (x) at the current iterate x k in order to generate the new iterate x k+1 . In our setting, however, the gradient information is known only locally at sensor i. Thus, the standard gradient descent method is not adequate.
To deal with the distributed nature of the sensor network information, we consider the incremental gradient method to minimize f (x), without the sensors explicitly sharing the functions f i (x) (see, [7], [13] and the references therein). In this algorithm, the iterates are generated according to The key difference between the standard gradient and incremental gradient method is that the standard gradient method generates iterates by using the gradient information of all functions f i (x) at the same (current) estimate x k , while the incremental method generates iterates through a cycle of intermittent adjustments z i−1,k+1 using only one function at a time, i.e., the gradient ∇f i (z i−1,k+1 ), so that all functions f i are processed within a cycle (see Fig. 1 for an illustration). The convergence of the incremental gradient method has been studied in [13], [7], [14] under different assumptions on the functions f i (x) and the step-size rules.

C. Recursive prediction error algorithm
Here, we discuss the standard recursive prediction error algorithm (RPE) for a parameter estimation problem (see [1]). To avoid confusion with the notation in the rest of the paper, we suppress the subscript i and consider the problem of estimating x from observations of a random process {R(k; x)} with the following dynamics: The RPE algorithm generates estimates of x by applying suitable approximations to the iterates generated by the gradient descent method as employed to solve an appropriate optimization problem. In particular, on the set X, the true parameter value x * minimizes the following function: When the standard gradient descent method is used to minimize f, the iteratesx k+1 are given bŷ The RPE algorithm obtains a sequence {x k } by using two approximations to the sequence {x k }.
The gradient of f (x) is not available in (10), but instead the sequence {r(k)} is available. The first approximation is a least mean-square (LMS) type approximation replacing the actual gradient ∇f (x k ) with an empirical gradient. Let us denote the iterate sequence corresponding to this approximation by The gradient ∇g k+1 (x; r k ), can be obtained from (5) and the extended representation of η (ℓ) k+1 (x; r k ) in (6). The problem is that even with this approximation the sequence {x k } cannot be obtained recursively.
Observe that to exactly evaluate g k+1 (x k ; r k ) and ∇g k+1 (x k ; r k ) one would need the entire vector r k .
To accommodate the recursive computations, we use another approximation Changing notation to reflect the approximations and re-ordering the equations, the resulting RPE algorithm can be written as follows, for ℓ = 1, . . . , d, The algorithm is initialized with values for ψ 1 , χ 1 and x 0 . Observe that to update x k the algorithm requires only r(k + 1), χ k+1 and ψ k+1 , and therefore, it is recursive. In summary, the iterates of the RPE algorithm are obtained from the standard gradient descent iterates with the following two approximations: 1) An LMS-like approximation for the gradient, and 2) An approximation to make the LMS approximations recursive.
The following theorem provides some sufficient conditions guaranteeing that the iterates generated by the RPE algorithm asymptotically converge to a minimum of f (x). The theorem is based on the results from [1].
Theorem 1: Let the following conditions hold.

1)
The set X is a closed and convex set containing x * . Furthermore, the system in (8) is stable, observable and controllable for all x ∈ X.
2) The matrices F (x) and G(x) are twice differentiable for all x ∈ X.
3) The fourth moments of V (k) are bounded. The second moments of W (k; x * ) are bounded.
Moreover, let the step-size α k be such that kα k → µ for some positive scalar µ. Then, the iterates x k generated by the RPE in (12) converge to a local minimum of f (x) in (9) over the set X, with probability 1.
Theorem 1 follows from Theorem 4.3 on page 182 and the discussions in pages 172 and 184 of [1].
The conditions for convergence of the algorithm are extremely weak. Note that the algorithm guarantees convergence only to a local minima and not necessarily to the global minimum when the function f (x) is convex this implies convergence to a global minimum.

VI. INCREMENTAL RECURSIVE PREDICTION ERROR ALGORITHM
As discussed in Section IV-A, when there are multiple sensors, the true parameter x * minimizes We combine the incremental gradient algorithm in (7) with the RPE algorithm in (12) to develop an incremental recursive prediction error (IRPE) algorithm. The main idea of the IRPE is to use an RPE like approximation for the gradient term in the incremental gradient algorithm (7). Formally, the iterates are generated according to the following relations for i ∈ I, and ℓ = 1, . . . , d, The initial values for the recursion are fixed at x 0 = x s , ψ i,1 = ψ i,s and χ i,s . To see that the algorithm has a distributed and recursive implementation assume sensor i − 1 communicates z i−1,k+1 to sensor i in slot k + 1. Sensor i then uses 4 r i (k + 1) to updates the iterate z i−1,k+1 to generate z i,k+1 . This is then passed to the next sensor in the cycle. Observe that in updating z i−1,k+1 , sensor i requires i,k+1 and ψ i,k+1 , which were calculated by sensor i in the previous time slot. Thus, the algorithm is recursive and distributed. Furthermore, note that sensor i only needs to know its own system

A. Convergence result
The iterates generated by the IRPE method are three approximations away from the iterates generated by the standard gradient descent method. The first approximation is in going from the standard gradient algorithm to the incremental gradient algorithm, the second is in approximating the gradient of the function with an LMS-like empirical gradient and the third is in calculating the empirical gradient recursively. Therefore, it is not clear if the iterates will converge to x * . We next state a theorem that provides sufficient conditions for the convergence of the iterates generated by the IRPE algorithm.
Theorem 2: For all i ∈ I, let the following conditions hold 1) The set X is a closed and convex set containing x * . Furthermore, the system in (8) is stable, observable and controllable for all x ∈ X.
2) The matrices F i (x) and G i (x) are twice differentiable for all x ∈ X.
3) The fourth moments of V i (k) are bounded. The second moments of W i (k; x * ) are bounded.
Moreover, let the step-size α k be such that kα k → µ for some positive scalar µ. Then, the iterates x k generated by the IRPE algorithm in (14)- (19) converge to a local minimum of f (x) in (13) over the set X, with probability 1.
Note that the result implies that for each i the iterates z i,k+1 converge to the same local minimum.
Thus the algorithm is not necessarily consistent.
There is alternative way to interpret the IRPE algorithm. In particular, consider a centralized scheme where the sensors immediately communicate their measurements to a fusion center. Now, at the fusion center, the RPE algorithm can be used to estimate the parameter x. For the specific model of our interest, there is a hidden structure in the RPE algorithm that permits an incremental implementation. Thus, the IRPE algorithm can also be viewed as an incremental implementation of a centralized RPE algorithm.
Since this hidden structure in the RPE algorithm is not easily identified, we have not used this alternative 4 We are assuming that sensor i obtains its measurement before it receives the iterate. From an implementation perspective, each time slot can be divided into two parts. In the first part, the sensors make measurements and in the second part they process. approach to actually present the algorithm. However, we use this approach to prove the convergence result stated in Theorem 2. The proof is provided in Appendix A.

B. Communication requirements
Incremental algorithm can potentially require less communication than centralized schemes. In a centralized scheme, in every slot each sensor has to communicate its measurements to a fusion center that is O(1) meters away on average. Summed over the m sensors in the network, the total communication requirement in a centralized scheme is O(m) bit meters per slot. In the incremental scheme, each sensor needs to pass only the iterate to a neighbor which is O log m √ m meters away on average, as discussed in [6]. Therefore, the total communication required in the incremental processing scheme is O( √ m log m) bit meters per slot.

C. Centralized versus incremental: Tradeoff
In our analysis we do not use any information about the joint statistics of the random process {Θ i (k; x)} and {Θ j (k; x)}. When this is the case the performance of the IRPE is identical to the performance of the centralized RPE algorithm.
Suppose some information about the join statistics is available. This information cannot be used in a distributed system because at most one sensor's measurement is known at a single location at any time.
Thus, the joint distribution information is not useful to the RPE.
A centralized system, on the other hand, can potentially use the joint density information to obtain a cost function f (x) that generates estimates with better properties. As an example, suppose that Θ i (k; x) = Θ j (k; x) for all k ≥ 1 and i, j ∈ I, which corresponds to the case when all sensors sense a field with no spatial variation, synchronously at time mk. Define H, respectively V (k + 1), to be the block matrix obtained by stacking the matrices H 1 , . . . , H m , respectively vectors V 1 (k + 1), . . . , V m (k + 1). Then, the centralized measurements R(k; x) have the following evolution: The corresponding time-invariant Kalman predictor is given by Notice that the predictor g i,k+1 (x; r k ) for the (k + 1)-st measurement of sensor i is a function of the past measurements made by sensor j, j = i. Using this predictor we can define a cost function in a manner similar to (3). As we will see in Section VIII, the nature of the cost function may be significantly better in terms of the number of local minima and the RPE applied to the system in (20) may have a better performance.
To summarize, there is an implicit tradeoff when we use the IRPE. Potentially better estimates may be obtained by a centralized scheme when the joint statistics of the process Θ i (k; x) and Θ j (k; x) are available. This is indicated by our numerical results in Section VIII.

VII. EXTENSIONS
We next discuss some extensions to the IRPE algorithm.

A. Hybrid scheme
Let us consider an alternative network architecture where the network of m sensors, divided into m c clusters of approximately equal size, is deployed in a unit square. Each cluster has a cluster head that is a neighbor to all the sensors in the cluster. We can develop a hybrid algorithm that is centralized intra-cluster and distributed inter-cluster. Each cluster head collects all the measurements made by the sensors in the cluster, and then the cluster heads use the IRPE algorithm to estimate x without sharing their measurements.
Note that, as each sensor is in the neighborhood of its cluster head, it is still required to only communicate to a neighbor. The cluster heads might have to communicate over larger distances. The that is available about the joint statistics of the processes seen by the sensors in the cluster.

B. Distributed and recursive regression
In the problem we have studied, the actual sensor measurement sequence {r i (k)} is a sample path of the random process {R i (k; x)} for x = x * . While we did not assume to know the value of x * , we did assume that for some x ∈ X the actual system is correctly modeled by (1).
In practice, it is very difficult to obtain the correct description and often approximate models are used.
In the context of our problem, this means that (1) need not necessarily be the correct description of the actual system dynamics for any value of x ∈ X. The minimum of f (x) is now interpreted as the value of x that chooses the state-space system that best approximates the actual system among all the systems generated as x ranges over X.
Theorem 4.3 of [1] concludes that even in this case, under some weak regularity conditions on the actual measurement sequence (instead of Condition 3) the iterates generated by the RPE still converge to a local minimum of f (x). Since the IRPE was proved to be equivalent to the centralized RPE, the above statement extends to the IRPE algorithm also.

C. Other extensions
We have not included an explicit input in modeling the system. Much of the analysis immediately follows where there is a deterministic open-loop input {u i (k)} that drives the system in (1). Of course, {u i (k)} should be known to sensor i. Another immediate extension is to the case when the matrix H i and noise V i (k) are also be parametrized by x. Finally, we remark that rate of convergence results are available for the RPE through a central limit theorem and these can be extended to the IRPE.

VIII. APPLICATION
We next consider a gas-leak problem to illustrate the concepts developed in the paper. We assume that a wireless sensor network is deployed inside a warehouse where gas tanks are stored. The network objective is to localize a leak, when one occurs. We use a two-dimensional model for this scenario, which is appropriate when the gas is significantly heavier than air. In any case, the extension to three dimension is immediate. We also remark that we have used the gas leak problem as only a representative example; the analysis is more generally applicable to heat and other diffusing sources.
Leak model: We asume that the leak occurs at time t = 0 and that the network detects the leak immediately. We model the leak as a point source at x = (x 1 , x 2 ). Each sensor sampling has a duration of 1 time unit. The leak intensity is modeled as a piece-wise constant function, i.e., the leak intensity is equal to I k during the time interval [k − 1, k) for k ≥ 1. Across sampling intervals, the leak intensity values vary according to the following Markov process: Here, ρ is a known scalar and {S(k)} is a sequence of i.i.d. Gaussian random variables with zero mean and variance σ 2 s . Thus, the intensity evolves in time as follows: where rect(t) is the rectangular function taking value 1 in the interval [0, 1] and zero elsewhere.
Medium model: We model the warehouse as a rectangular region with known dimensions l 1 × l 2 .
Without loss of generality, we let the warehouse to be the region D = [0, l 1 ] × [0, l 2 ], and we denote the boundary of the warehouse by ∂D.
The medium is characterized by the diffusion coefficient of the gas, boundary conditions and initial conditions. We use C(y, t; x) to denote the concentration at a point y at time t when the source is at x.
We make the following assumptions.
1) The diffusion coefficient of the gas is the same everywhere in the warehouse. We use ν to denote this value.
Observation model: Let s i be the location of the i-th sensor. We assume that all sensors sense at the beginning of each time slot, i.e., at time k. Then where N i (k) is a zero mean i.i.d. measurement noise with known variance σ 2 n . Transport model: We assume that the transport of the gas in the warehouse obeys the diffusion equation.
Here, ∇ 2 is the Laplacian differential operator andδ is the Dirac delta function.

A. Problem statement and related literature
The medium characteristics are completely known, i.e., l 1 , l 2 , and ν are known. The sensors' sampling duration and the measurement noise variance σ 2 n are also known. Further, the variance σ 2 s of S(k) is known. The problem is to determine the location of the point source x from the sensor measurements in a distributed and recursive manner. To solve the above problem we first show that, as a consequence of the assumptions that have been made, the sensor measurements follow a state-space model. We then use the IRPE algorithm developed in the previous sections to solve the problem.
We next compare and contrast the models described above with the models used in literature. The point source model is a common model for diffusing sources and has been extensively used in localization studies [15]- [18]. The random time-varying source intensity model is more realistic compared to the constant intensity [15], [17], [18] and instantaneous intensity models that are usually studied. Localization of sources with time-varying intensity have been studied in a centralized and non-recursive setting in [19], [20]. These studies consider a deterministic evolution of the leak intensity and use a continuous observation model. We are not aware of any paper that models the time-varying intensity as a random process. Most papers study that case when the medium is infinite or semi-infinite since the diffusion equation has a closed form solution in that case [15], [17]. The medium model assumed in this paper is more general. We also remark that we can extend the results to non-rectangular geometries by using the Galerkin approximation [21].
While centralized recursive source localization has recieved much interest [15], [17], [20], [22] there are very few papers that discuss a distributed solution. A recursive and distributed solution to the problem in a Bayesian setting is discussed in [16]. A related paper is [23] that deals with the problem of estimating the diffusion coefficient in a distributed and recursive manner. We are not aware of any prior work that solves the source localization problem using a distributed and recursive approach in a non-Bayesian setting.

B. Approach
We show in Appendix B that by using Green's technique to solve differential equations it is possible to obtain a state-space description for each sensor's observation process. We can then use the IRPE algorithm to estimate the iterates in a distributed and recursive manner.
The initial intensity value is taken to be 100, ρ is fixed at 0.99 and the variance of S(k) is fixed at 10.
A network of 27 sensors is deployed. To ensure complete coverage of the sensing area, we first placed 9 sensors on a grid and then randomly deployed 2 sensors in the immediate neighborhood of each of the 9 sensors. The network is shown in Fig. 2.  Centralized RPE with information IRPE with cluster size 3 Actual value IRPE with cluster size 1 processing through the IRPE algorithm this information is not used.
3) Centralized RPE: All sensors immediately communicate their measurements to a fusion center. In this case the information is completely used. The fusion center runs 1000 iterations.
The results are plotted in Figs. 3 and 4. As expected the centralized RPE performs the best. However, what is interesting to note is that the standard IRPE does not converge to the correct solution but is caught in a local minimum. We also observed this in other simulation runs. However, when the sensors are clustered the iterates converge to the correct location.

IX. CONCLUSIONS
Linear state-space models arise naturally, or as linear approximations to non-linear state-space models, in many applications. In an inference setting where the aim is to estimate a quantity of interest, it is quite natural for the state-space models to be parametrized by the unknown quantity of interest. In a control system setting where the aim is to control the process, it is common to use such incomplete state-space models as 'grey-box' descriptions of the system that is to be controlled. Thus, the problem addressed in this paper is important in both of these settings.
In Section VII we could only give a qualitative description of the tradeoff between centralized and distributed schemes. It is therefore of interest to find good bounds on the performance of the IRPE and RPE schemes that can be used to quantify the loss in performance. Also, to truly understand the performance of the algorithm in practical settings, we need to obtain convergence results when there are communication errors. Further, we have considered a simple class of networks where the topology is fixed. It is important to obtain an algorithm that is similar to the IRPE for networks with a random and time-varying topologies. Finally, as mentioned in Section VII, the result and analysis extend easily to the case where there is an open-loop input to the system. An important extension is to obtain similar convergence results for some common classes of closed-loop inputs using the techniques discussed in Appendix 7.A of [1].

A. Proof of Theorem 2
We take the following approach to prove Theorem 2. First, we consider a centralized system where the sensors immediately communicate all their measurements to a fusion center. For this system we use the RPE algorithm to generate a sequence of iterates, and then we show that this iterate sequence is identical to the iterate sequence generated by the IRPE algorithm. We complete the proof by proving that the iterates converge to a local minimum of f (x).
In what follows, we extensively use the notion of block vectors and matrices. For positive integers a and b, let M a×b be the vector space of all real matrices of dimensions a × b. A block vector in M a×b is a vector whose elements are from M a×b . The length of a block vector is the number of block elements.
In a similar manner, block matrices in M a×b are matrices where each element is itself a matrix from M a×b . While writing block matrices we will allow for a slight abuse of notation and use 0 and I to denote the zero and identity matrices, respectively. Their dimensions can be unambiguously fixed from the dimensions of the other blocks in the block matrix. We will use U a b , b ≤ m, to denote the unit block vector in M a×a of length m, with the b-th block equal to the identity matrix in M a×a and all the other blocks equal to the zero matrix in M a×a .
We allow i, j to take values in the set I = {1, . . . , m}. We define δ[·] as the Kronecker delta. Recall that the dimension of the matrices Θ i (k; x) is q, the dimension of the measurement r i (k) is p, and the dimension of the parameter vector x is d. Also, recall that for any random process {Y (k; x)} that is parametrized by x, y(k) denotes the sample path of Y (k; x * ). Denote the sequence communicated by a sensor i by {r i (mk + j)}, with Next, denote the observation sequence at the fusion center by {r(mk + j)}, wherẽ We now consider the problem of estimating x from observation sequence {r(mk + j)}. We show that the random processR(mk + j; x) can be represented as the output vector of a suitably defined statespace system. For this, we first use the relations in (1) to obtain the equations describing the evolution {R i (mk + j; x)}. Note that from (25), we havē LetD i (x) be the following m × m block matrix in M q×q : Observe thatD Also, defineH i = H i (U q 1 ) T , and note that The following is an illustration forΘ i (mk + j; x) with i = 3 and j = 2, 3, 4: Claim 1: For all n ≥ 0, we havē Proof: Let n = mk + j. Substitute forΘ i (n; x) andW i (n; x) from (31) in the right hand side (RHS) of (33). For i < j, from (29) we obtain RHS of (33 For i = j, using (29) and (1), we obtain Finally, when i < j, from (29) we have thus, completing the proof of the relation in (33).
We next prove the relation in (34). At first, we consider the case when j = m and show that the following relation holds: Let i = j + 1, and note that from (32) we haveV i (n + 1) = 0. Furthermore, from the definition of Using the expression in (30), we can immediately verify thatH iΘi (mk + (j + 1); x) = 0. Therefore, we haveR i (n + 1; x) = 0 for i = j + 1.
When j = m, by using arguments similar to that of the preceding case j = m, we can show that thus, completing the proof. Now, by combining the equations in (33)-(34) for i ∈ I, we provide evolution equations for {R(n; x)}. DefineF Using the relations in (33) and (34), we can writẽ R(n + 1; x) =HΘ i (n + 1; x) +Ṽ (n + 1).
To apply the RPE we need to determine a predictor family forR(n; x * ) that is parametrized by x and is asymptotically optimal at x = x * . We do this in the next section.
As we mentioned earlier,r(n) contans the same information asr(n) about the true value of x. Define for ℓ = 1, . . . , d, We define the iterates {x n } as follows: Here, α(n) = α k+1 for n = mk + j for j = 1, . . . , m − 1. Next, we assign the initial values for the recursion. Recall that the IRPE algorithm in (19) is initialized with the values ψ i,1 = ψ i,s , ξ i,s for all i and ℓ, and x 0 = x s . We letx 0 = x s , and i,s for all i and l. 4) Rest of the proof: Finally, here we show thatx n = z j,k+1 . Recall that ψ i,k and χ (ℓ) i,k are generated in the IRPE algorithm (14)- (19). Define for l = 1, . . . , d, We next establish the following lemma. Once we prove this lemma we can use induction to show that the iterates generated by the IRPE algorithm are the same as those generated by the centralized scheme.
When i < j + 1, we haveψ Equation (53) can be shown similarly for the case j = m. The proof of relation (54) is very similar to that of relation (53) and therefore, it is omitted.
Observe that the initial values for the RPE algorithm in (49) and the IRPE algorithm in (19)  By using Lemma 1 and the induction on k, we can conclude thatx mk+i = z i,k+1 for all k ≥ 1 and i ∈ I. To complete the proof, we only need to show that the sequence {x n } converges to a minimum of f (x), which is done in the following.  (13), over the set X w.p.1.
Proof: By the assumptions of Theorem 2, it follows that the conditions of Theorem 1 are satisfied.
Therefore, by Theorem 1, the iteratesx n converge to a local minimum over the set X of the following functionf (x) = lim N →∞ 1 N N n=1 E R (n; x * ) −g n (x,R n (x * ))

B. State-space Model for the Sensor Measurements
A standard technique to solve partial differential equations with boundary and initial conditions is to use Green's function. For the equation in (24) the solution can be written as follows Here, τ and z = (z 1 , z 2 ) are parameters of integration, and G f is the following Green's function where To get a convenient approximation we will use a sufficiently large, but fixed number of terms in the convergent infinite series in (57). We will let n i , i = 1, 2, vary from 1 ton i , where the integersn i > 0 are chosen large enough to provide a sufficiently good approximation. Therefore, C(y, t; x) ≃ cos n i πx i l i , Θ ′ n1,n2 (t; x) = A ′ (x, n 1 , n 2 ) t 0 I(τ )β t−τ n1,n2 dτ .
From (21), (60) and (23)  Hence, complete information is available about the joint statistics of Θ i (k; x) and Θ j (k; x). If the sensors do not sense synchronously at the beginning of the slot, the model is more complex. Nevertheless, we can still identify a state vector and write the sensor measurements as the output process of a state-space system.