Sandy Mahfouz
email: sandy.mahfouz@utt.fr

Farah Mourad-Chehade

Paul Honeine
email: paul.honeine@utt.fr

Joumana Farah
email: farah.mourad@utt.fr

Hichem Snoussi
email: hichem.snoussi@utt.fr

Gas Sources Parameters Estimation Using Machine Learning in WSNs

Keywords: gas diffusion, machine learning, one-class classification, ridge regression, source parameter estimation I

This paper introduces an original clusterized framework for the detection and estimation of the parameters of multiple gas sources in WSNs. The proposed method consists of defining a kernel-based detector that can detect gas releases within the network's clusters using concentration measures collected regularly from the network. Then, we define two kernel-based models that accurately estimate the gas release parameters, such as the sources locations and their release rates, using the collected concentrations.

estimates. Another method was proposed in [START_REF] Delmaire | Joint estimation decision methods for source localization and restoration in parametric convolution processes. application to accidental pollutant release[END_REF], where the problem of pollutant source localization and flow estimation is addressed in a one-dimensional context, using a single remote sensor. The pollutant is assumed to be generated by one out of several possible sources, and the task is viewed as a conditional deconvolution which requires a priori knowledge. In the end, a joint estimation decision is derived in a Bayesian framework.

In terms of computational concept, algorithms in WSNs can be divided into centralized, distributed, and clusterized algorithms [START_REF] Nakamura | Information fusion for wireless sensor networks: Methods, models, and classifications[END_REF]. The centralized algorithms require the transmission of all measurements to a fusion center (e.g., a sink node) for processing. Such strategies often result in prohibitive wasteful energy and bandwidth consumptions and can thereby reduce the lifetime and utility of the network. In this paper, we propose a new clusterized framework for the detection and estimation of multiple gas sources in wireless sensor networks. Sensors are uniformly or randomly deployed in the region of interest and collect concentration measurements at regular and short sampling intervals. The region is divided into clusters, each having its own cluster head. Consequently, information processing is done locally. Compared to a centralized strategy, this clusterized approach is more robust to failures, since several cluster heads are engaged in the detection phase. It is also less energy consuming and thus more adapted to WSNs limitations.

The proposed framework consists of two phases: the detection phase and the estimation phase. In the first one, a kernel-based detector is defined using the Support Vector Data Description (SVDD) [START_REF] Tax | Support vector data description[END_REF]. The detector is used by each cluster head in order to identify any gas release in the specified cluster. When an anomaly is spotted, the concentration vector that first triggered the alert is treated to estimate the gas release parameters, such as the source's location and the release rate. Note that not all the concentrations will be processed in the estimation phase, since some are irrelevant for parameters estimation. Details about the selection of the useful information will be given in the sequel. A first estimate of the parameters is obtained using a non linear model, which is also defined within the framework of kernel methods, proved to be successful for solving nonlinear regression problems [START_REF] Mahfouz | Kernel-based localization using fingerprinting in wireless sensor networks[END_REF], [START_REF] Mahfouz | Ridge regression and kalman filtering for target tracking in wireless sensor networks[END_REF], [START_REF]Target tracking using machine learning and kalman filter in wireless sensor networks[END_REF]. Then, part of the estimated parameters are processed as internal feedback, along with the measured concentrations, in order to provide a more accurate estimate of the source location. Simulations show that the proposed estimation method yields accurate results in the case of a single source, as well as in the case of multiple sources.

The rest of the paper is organized as follows. Section II introduces the considered advection-diffusion model. Section III provides a thorough description of the detection phase in the proposed framework. In Section IV, the estimation phase is described, along with details about the definition of the two kernel-based regression models. Next, in Section VI, the effectiveness of the method is discussed for different scenarios, and a comparison to other methods is provided. Finally, Section VII concludes the paper.

II. THE ADVECTION-DIFFUSION MODEL

This section describes the advection-diffusion model considered in this paper for the generation of the concentrations. Note that this model is a generalized version of the one developed in [START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF]. Now consider an instantaneous gas release of Q kg, assumed to occur at time t 0 at location (x 0 , y 0 , z 0 ). A wind with mean velocity U = (u x , u y , 0) spreads the gas particle in the region. The mass concentration C of the released agent, at an arbitrary location (x, y, z) and at time t, is governed by the equation of mass conservation given by the following:

∂C ∂t = -∇q, (1)
where ∇ is the gradient operator, and q is the pollutant mass flux per unit area. The flux q is given by:

q = C U -   K x 0 0 0 K y 0 0 0 K z   ⊗ ∇C, (2)
where C U is the mean mass advection by the wind, ⊗ is the tensor product, and K x , K y , K z are eddy diffusivities in the X, Y and Z directions respectively. Equation (2) can then be written as follows:

q = C u x -K x ∂C ∂X , C u y -K y ∂C ∂Y , -K z ∂C ∂Y . (3)
By substituting (3) into (1), we get an equation that can be solved subject to two boundary conditions [START_REF] Christopoulos | Multi robot trajectory generation for single source explosion parameter estimation[END_REF], [START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF]. The first condition results from the fact that the concentration is zero at infinity in all spatial directions, and the second condition is that the gas is not absorbed by the ground. To simplify the model, the velocity of the wind U , as well as the eddy diffusivities K x , K y , K z , are assumed to be constant. Following these assumptions, we get the following solution:

C(x, y,z, t) = Q 8 π 3 2 (K x K y K z ) 1 2 ∆t 3 2 × exp - (∆x -u x ∆t) 2 4K x ∆t - (∆y -u y ∆t) 2 4K y ∆t × exp - ∆z 2 4K z ∆t + exp - ∆z ′ 2 4K z ∆t , (4)
where ∆x = xx 0 , ∆y = yy 0 , ∆z = zz 0 , ∆z ′ = z + z 0 and ∆t = tt 0 .

In the following, for the sake of simplicity, we assume that all measurements are taken at ground level. We also suppose that the gas release occurs at ground level, which means that the gas source is at location (x 0 , y 0 , z 0 ) = (x 0 , y 0 , 0). Therefore, using (4), the concentration measured at a location (x, y, 0) is given by the following:

C(x,y, 0, t) = Q 4 π 3 2 (K x K y K z ) 1 2 ∆t 3 2 × exp - (∆x -u x ∆t) 2 4K x ∆t - (∆y -u y ∆t) 2 4K y ∆t . (5)
Note that the wind velocity U could be provided by an anemometer, and is therefore treated as a known constant [START_REF] Christopoulos | Multi robot trajectory generation for single source explosion parameter estimation[END_REF]. The other parameters of the model are unknown and need to be estimated when a gas explosion occurs over the area under scrutiny. Hence, having gas concentrations measured over a certain area using a WSN, the aim of this paper is to detect the gas explosion when it occurs. Then, the objective is to estimate the location (x 0 , y 0 , 0) of the explosion, the gas release mass Q, and the eddy diffusivities K x , K y and K z , with K x and K y assumed to be equal [START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF].

III. THE DETECTION PHASE

As we already stated, the proposed method uses a WSN to continuously measure gas concentrations over a certain area Ω. Then, computations consist of two phases: the detection phase and the estimation of the source's parameters phase. In this section, we first give a brief description of the network's setup and of the detection phase. Then, we define a classifier that is capable of optimally separating the data into normal and abnormal data, where the abnormality denotes the occurrence of a gas diffusion.

A. Description of the detection phase

Consider N sensors deployed in the area Ω to be monitored, at fixed locations (x n , y n , z n ) = (x n , y n , 0), n ∈ {1, . . . , N }. These sensors measure, at each time t, the gas concentrations at their locations. Let (C(x 1 , y 1 , z 1 , t) . . . C(x N , y N , z N , t)) ⊤ be the N × 1 vector of the gas concentrations measured at time t by the N sensors. These concentrations are assumed to follow the advection-diffusion model, detailed in Section II. In order to thoroughly monitor the area of interest, the concentrations are measured at regular and short sampling intervals by the deployed sensors.

Assume now that an instantaneous gas release of Q kg occurs at time t 0 at location (x 0 , y 0 , z 0 ) = (x 0 , y 0 , 0), all these parameters being unknown. The gas particles are then spread by a wind with a mean known velocity U = (u x , u y , 0). Since a gas is released in the area, an alarm should be triggered based on the vector of measured concentrations. To this end, we need to develop a detector, capable of determining whether the concentrations are normal or not.

We propose to partition the area into Z distinct clusters, in order to make the computations less complex and the method more robust to transmission impairments and network failures. Each cluster is managed by a cluster head, that is a smart central processing unit (CPU), responsible of handling (gathering and synchronizing) data, performing calculations, and exchanging information with the sensors. Note that this device could also be one of the sensors of the network, and it can be placed anywhere in its specified cluster. In addition, all cluster heads can communicate with each other, and clusters could have any shape or dimension. Without loss of generality, we consider here that clusters are rectangular, with the cluster heads located at their centers, as illustrated in Fig. 1. Now that the network has been configured, each cluster head receives the measured concentrations from the sensors in its cluster at every time t. Let n z denote the number of sensors in cluster z, z ∈ {1, ..., Z}, and let C (z) (t) denote the concentration vector of size n z ×1 collected by the cluster head of cluster z at time t. The received data are then processed instantly to detect whether a gas diffusion has occurred or not.

To achieve this, we aim at finding a detector that can detect an anomaly in the monitored region; in other words, it can detect whether or not there is one or more gas diffusions. A one-class classifier is able to do this task. Indeed, the Support Vector Data Description (SVDD) has been introduced to address the problem of anomaly detection [START_REF] Tax | Data domain description using support vectors[END_REF], [START_REF] Tax | Support vector data description[END_REF].

B. Definition of the detector using SVDD

The objective here is to define the class boundary, given data that are originated from a single class, i.e., the normal data, but are possibly contaminated with a small number of outliers, i.e., the abnormal data. For simplicity, the number of sensors in each cluster is assumed to be the same, so that the same defined detector could be used for all the clusters, as will be explained in the experimental results. As we already explained, a one-class classifier can be used for the detection, such as the Support Vector Data Description (SVDD) [START_REF] Schölkopf | Introduction to support vector learning[END_REF], [START_REF] Schölkopf | A generalized representer theorem[END_REF], [START_REF] Tax | Support vector data description[END_REF]. It essentially fits the smallest possible sphere around the given normal data, allowing some samples to be excluded as outliers. Therefore, a spherically shaped decision boundary with minimum radius is computed to enclose most of the training data. Data lying outside this decision boundary are considered as abnormal, i.e., outliers.

Let φ denote a function that maps the data from the input space IR nz , where the concentrations lie, into a higher dimensional feature space H. Consider a reproducing kernel κ : IR nz ×IR nz → IR, with H its reproducing kernel Hilbert space (RKHS) with inner product •, • H . We then have

κ(C i , C j ) = φ(C i ), φ(C j ) . Now consider a training set C i , i ∈ {1, .
.., N det }, N det being the size of the training set for the detection phase. The SVDD consists of estimating the hypersphere with minimum radius that encloses all data φ(C i ) in the feature space H. Let a be the center of the hypersphere, and R > 0 its radius. To allow a better description of data, we allow the presence of outliers in the training set by introducing the slack variables ξ i ≥ 0. Consequently, abnormal concentration vectors are also used in the training set. Mathematically, the values for the slack variables are obtained by minimizing a cost function that balances the volume of the hypersphere against the penalty associated with outliers. Note that minimizing the hypersphere's volume is equivalent to minimizing R 2 . Therefore, we get the following constrained optimization problem:

min a,R,ξi R 2 + 1 νN det Ndet i=1 ξ i (6) subject to φ(C i ) -a 2 H ≤ R 2 + ξ i and ξ i ≥ 0 ∀i = 1, ..., N det . (7)
The quantity ν is a predefined parameter that regulates the trade-off between the volume of the hypersphere and the number of outliers.

Let L denote the Lagrangian of the above constrained optimization problem. By taking the partial derivatives of L with respect to R, a and ξ i , we get the following relations:

Ndet i=1 α i = 1, a = Ndet i=1 α i φ(x i ), and 0 ≤ α i ≤ 1 νN det ,
where α i are the Lagrangian multipliers. Incorporating these relations into the Lagrangian L gives us the following objective functional to be maximized with respect to α i :

L = Ndet i=1 α i κ(C i , C i ) - Ndet i=1 Ndet j=1 α i α j κ(C i , C j ), (8)
subject to 0 ≤ α i ≤ 1 νNdet
. This is a quadratic programming problem, whose solution is found using an off-the-shelf optimization technique. For instance, one can use the Matlab function quadprog to solve such a problem, and consequently compute the α i .

As for the radius of the optimal hypersphere, it is, by definition, the distance from the center a to any sample φ(C k ) on the boundary in the feature space H. Therefore, the radius is given by the following:

R 2 =κ(C k , C k ) -2 Ndet i=1 α i κ(C k , C i ) + Ndet i=1 Ndet j=1 α i α j κ(C i , C j ).
Sensors record the concentrations at regular short intervals T . Therefore, each T seconds, a new concentration vector C (z) (t d ) is taken from each cluster z, with

t d = d × T, d ∈ IN.
The decision rule for each new concentration vector C (z) (t d ) is obtained by evaluating the distance between the center a and the mapping φ(C (z) (t d )) in the feature space, given by φ(C (z) (t d ))a 2 H . Note that to compute this distance, one does not need the exact expressions of φ and a. Instead, by developing the distance expression, one gets scalar products of φ functions, that could be replaced by kernels. The new concentration

C (z) (t d ) is considered as normal if the distance calculated is smaller than the radius, i.e., φ(C (z) (t d )) -a 2 H ≤ R 2 .
Otherwise, an alert is triggered in the considered cluster. According to [14, Section 7.1], we can provide an upperbound on the false alert rate, that is the probability that a normal concentration vector is mistakenly classified as abnormal.

The proposed algorithm can detect and estimate more than one source. Therefore, more than one detector might detect a gas release at the same time. However, for simplicity, all the diffusions are assumed not to overlap. It is worth noting here that by dividing the whole area into several clusters, the computations are processed in parallel at all the cluster heads. Moreover, the size of the considered data in the classifier for training and detection is limited to n z , instead of taking all the sensors' concentrations at once. This way, the network could be as large as needed without increasing the computations' complexity.

IV. THE ESTIMATION PHASE

Now that our detector is defined, we are able to detect a gas diffusion in a specified cluster. Once an alert is triggered in a cluster z, we proceed to the processing of the concentration vector C (z) (t) recorded by the sensors in the specified cluster. In this section, we give a brief description of the estimation phase, and we emphasize on the choice of the concentrations to be used in the source's parameters estimation. The idea here consists of grouping the clusters in alert in a way to discern all the sources of explosion. One can take all the clusters in alert together; however, if the network is too large, with many sources of explosion, one gets high dimension data to process. An alternative way consists of grouping adjacent clusters in alert together; by adjacent, we mean the clusters that share a common boundary. Each cluster in alert but not adjacent to others will form a group on its own. To illustrate this, see Fig. 2, where six explosions are assumed to occur in the region of interest.

At the end of this step, several groups with different sizes are generated, each covering one or more explosion sources. The concentration vectors recorded by the sensors of a group are then communicated to a cluster head of the group, e.g., the one having the highest computation capabilities, and the concentrations are concatenated in a single group vector. Now, in order to detect the multiple gas releases per group, the proposed method uses the concentration group Fig. 2: Region on interest in the case of multiple gas releases-represents the cluster heads when there is no alert, when there is an alert, * represents the locations of explosions, and the red rectangle represents the defined groups. Grouping the clusters in alert according to their adjacencies leads here to four groups. * Fig. 3: Group topology-represents the cluster heads, • represents the deployed sensors, * represents the location of the explosion, and represents a local maximum. vectors to compute local maxima, leading ideally to the concentrations of the closest sensors to the gas sources. Let N s be the total number of local maxima detected at all groups. This means that N s sources have exploded at time t, or have exploded at different times, but are being detected and estimated simultaneously, since the maximal concentrations would be at the release source. Consider now the example of Fig. 3, where the source of explosion is near the boundary of two clusters. The alert is then triggered in both clusters, the gas being spread all around the source. Since the clusters are adjacent, they are grouped together in order to form one group. Then, the local maximum is estimated, leading to only one sensor, instead of two if the clusters were considered separately. This means that by grouping the clusters, one would identify the exact number of exploded sources, even if the gas is spread over several clusters.

Having detected the local concentrations maxima, we now need to process the information in each group in order to estimate all N s sources' parameters. It is important to mention that not all concentration information have the same importance in the algorithm. Consequently, only the information collected around the local maxima will be treated. This solution relies on the fact that the highest con-centrations are collected by the sensors around the source of explosion, and these concentration measures are the most relevant for parameters' estimation. If we reconsider Fig. 3, only the concentration information from the sensors around the red triangle denoting the maximum will be considered. Note that the size of the clusters and the number of sensors in the figure are only chosen for illustrative purposes.

In this paragraph, we propose to choose the concentrations to be treated based on the advection-diffusion model of Section II. In other words, we need to estimate the size of the small local zone around the local maximum as illustrated in Fig. 3. Let C thres be a threshold under which the concentrations are considered really small, which means that they are measured by sensors far from the source. Depending on the considered application and the type of monitoring needed, one can determine an approximation of the maximum release mass Q max . Using C thres and Q max , we can find the maximal distance between a source and a sensor. Then, this distance is used to determine the boundaries of the local zone, whose center is the local maximum and the estimated maximal distance is the distance separating the local maximum from the boundary of the zone. By taking the logarithm of (5), we get the following:

log C thres = log Q max 4 π 3 2 (K x K y K z ) 1 2 - 3 2 log ∆t - ∆x 2 + u 2 ∆t 2 -2u∆t∆x 4K x ∆t + ∆y 2 4K y ∆t
For simplicity, assume now that the local zone is circularly shaped, which leads to ∆x = ∆y. In addition, having K x = K y yields the following:

log C thres = log Q max 4 π 3 2 (K 2 x K z ) 1 2 - 3 2 log ∆t - 2∆x 2 4K x ∆t - u 2 ∆t 4K x + 2u∆x 4K x .
Finally, a second degree polynomial in terms of ∆x is obtained as follows:

a∆x 2 + b∆x + c = 0, (9)
where

           a = -2 4Kx∆t , b = 2u 4Kx , c = log Qmax 4 π 3 2 (K 2
x Kz)

1 2 -3 2 log ∆t -u 2 ∆t
4Kx -log C thres . The local zone dimensions ∆x = ∆y are then obtained by resolving the quadratic equation ( 9), using Q max , C thres , and initial approximative values of K x and K z . Then, only the concentrations measured by the sensors in that zone will be considered for the source parameters estimation; in other words, the considered concentrations are the ones measured by the sensors at a maximal distance ∆x from the local maximum. For simplicity, we assume that we have the same Q max and C thres for all the groups, thus for all N s local maxima. Therefore, the same local zone size will be considered for all N s explosions. Nevertheless, depending on the type of application, one can define different sizes for the local zones around the N s maxima. Now let c (s) (t), s ∈ {1, ..., N s }, denote the vector of useful concentrations collected at time t around each local maximum s, i.e., around each source s. For the sake of clarity, we drop the time t in the following, since only the concentrations around the detection time are used. The first objective is to define a model ψ A that takes as input the concentration vector c (s) and yields as output the vector of parameters θ (s) of source s, which includes the gas release mass Q (s) , the source location (x (s) 0 , y (s) 0 , z (s) 0 ) and the diffusivity constants. Then, the next objective is to define a second model ψ B that provides an enhanced estimation of the source location. Kernel methods [START_REF] Hofmann | Kernel methods in machine learning[END_REF] provide an elegant framework to define both models, as it will be shown in the following section. Note that estimations are handled by the cluster heads that are the closest to the local maxima sensors.

V. DEFINITION OF THE MODELS ψ A AND ψ B

Let V be the number of sensors in the local zones. As we mentioned earlier, only the concentrations around the local maxima will be considered in the estimation phase. Therefore, the vector c (s) of size V × 1 is the vector that will be processed to find the parameters of source s. In the following paragraph, we define the model ψ A , that takes as input the concentration vector c (s) and yields as output the source's parameters θ (s) , that are the gas release mass Q (s) , the source location (x (s) 0 , y (s) 0 , 0) and the diffusivity constants. Then, we define a new model ψ B that provides an enhanced location estimation of the source, using a part of the estimated source's parameters and the measured concentration vector c (s) .

A. First estimation

In this subsection, we define a model ψ A : IR V → IR 5 , that associates to each concentration vector the corresponding source's parameters. We propose to find the model ψ A by solving a nonlinear regression problem. The main benefit of such an approach is that no prior knowledge of the model is needed. Kernel methods [START_REF] Hofmann | Kernel methods in machine learning[END_REF] have been remarkably successful for solving nonlinear regression problems. More specifically, we consider the kernel ridge regression [START_REF] Saunders | Ridge regression learning algorithm in dual variables[END_REF] to determine ψ A , by combining five separate optimization problems, one for each component of the output vector.

Consider the following training set (c ℓ , θ ℓ ), ℓ ∈ {1, . . . , N reg }, where N reg is the size of the training set used for the regression phase and θ ℓ = (Q ℓ K x K z x 0ℓ y 0ℓ ). The vector c ℓ yields the concentrations recorded by V sensors uniformly distributed in a zone of the same size as the small local zone introduced previously, around the time of a gas release of parameters θ ℓ , with (x 0ℓ , y 0ℓ , 0) being the source's location and Q ℓ the gas release mass. As already mentioned, the eddy diffusivities K x and K y , with K x = K y , are assumed constant, since they depend on the atmospheric conditions and the type of chemical agent, considered constant as well. Note here that the advantage of using a local small zone to select the relevant concentrations allows us to save computations and reduce the training complexity, since V sensors are now considered instead of N . In the following, let Θ = (θ ⊤ 1 . . . θ ⊤ Nreg ) ⊤ . The matrix Θ is then of size N reg ×5, having Θ ℓ,i for the (ℓ, i)-th entry, with ℓ ∈ {1, . . . , N reg } and i ∈ {1, . . . , 5}. We also denote θ ℓ by Θ ℓ, * , and the i-th column of Θ by Θ * ,i .

Let ψ A = (ψ A1 . . . ψ A5 ), where ψ Ai , i ∈ {1, . . . , 5} estimates Θ ℓ,i , the i-th component of the vector Θ ℓ, * , for an input c ℓ . Each function ψ Ai is then determined by minimizing the mean quadratic error between the model's outputs ψ Ai (c ℓ ) and the desired outputs Θ ℓ,i : [START_REF]Target tracking using machine learning and kalman filter in wireless sensor networks[END_REF] where η A is a regularization parameter that controls the tradeoff between the training error and the complexity of the solution. According to the representer theorem [START_REF] Hofmann | Kernel methods in machine learning[END_REF], the optimal function ψ Ai can be written as follows:

min ψAi∈HA 1 N reg Nreg ℓ=1 (Θ ℓ,i -ψ Ai (c ℓ )) 2 + η A ψ Ai 2 HA ,
ψ Ai (•) = Nreg ℓ=1 β ℓ,i κ A (c ℓ , •), (11)
where " • " is the function's input, κ A : IR V × IR V → IR is a reproducing kernel, and β ℓ,i , ℓ ∈ {1, . . . , N reg }, are parameters to be determined. We now denote by β the N reg × 5 matrix whose (ℓ, i)-th entry is β ℓ,i . The vector β * ,i denotes then its i-th column, and β ℓ, * its ℓ-th row. By injecting [START_REF] Tax | Data domain description using support vectors[END_REF] in [START_REF]Target tracking using machine learning and kalman filter in wireless sensor networks[END_REF], the dual optimization problem in terms of β * ,i is obtained, whose solution is given by taking the derivative of the corresponding cost function with respect to β * ,i and setting it to zero. We then obtain the following:

β * ,i = (K A + η A N reg I) -1 Θ * ,i ,
where I is the N reg × N reg identity matrix, and K A is the N reg × N reg Gram matrix whose (v, w)-th entry is κ A (c v , c w ), for v, w ∈ {1, . . . , N reg }. For an appropriate value of the regularization parameter η A , the matrix between parenthesis is always non-singular. Since the same matrix (K A + η A N reg I) needs to be inverted in order to estimate each source's parameter, all five estimations can be collected into a single matrix inversion problem, thus reducing the computational complexity, as follows:

β = (K A + η A N reg I) -1 Θ. ( 12
)
Finally, using equation [START_REF] Tax | Data domain description using support vectors[END_REF] and the definition of the vector of functions ψ A (•), we can write ψ A as follows:

ψ A (•) = Nreg ℓ=1 β ℓ, * κ A (C ℓ , •). ( 13
)
Having detected a local maximum in the region, with its relevant concentration vector c (s) , one could then obtain a first estimate of the source s parameters as follows:

( Q Kx Kz x0 ŷ0 ) = ψ A (c (s) ). ( 14
)

B. Enhancement of the estimates

In this subsection, we introduce the second model ψ B , that takes as input the measured concentrations, the estimated mass release and the estimated eddy diffusivities, and gives as output an enhanced estimate of the location of the source. Using the already estimated information as internal feedback is expected to provide an improved accuracy of the first location estimates of the sources.

The definition of the model ψ B is also done using the ridge regression [START_REF] Saunders | Ridge regression learning algorithm in dual variables[END_REF]. However, the training set is defined differently, such that the training input is given by

W ℓ = (c ⊤ ℓ Q ℓ K x K z ) ⊤ and the training output is the source's location (x 0ℓ y 0ℓ ), ℓ ∈ {1, . . . , N reg }. Let ψ B = (ψ B1 ψ B2 )
, where ψ B1 and ψ B2 estimate x 0ℓ and y 0ℓ respectively. Let X 0 denote the N reg × 1 vector whose ℓ-th entry is given by x 0ℓ , and Y 0 denote the N reg × 1 vector whose ℓ-th entry is given by y 0ℓ .

Here as well, the models ψ B1 and ψ B2 are obtained by minimizing the mean quadratic errors between the estimated outputs and the desired ones, as follows:

min ψB1∈HB 1 Nreg Nreg ℓ=1 (x 0ℓ -ψ B1 (W ℓ )) 2 + η B ψ B1 2 HB min ψB2∈HB 1 Nreg Nreg ℓ=1 (y 0ℓ -ψ B2 (W ℓ )) 2 + η B ψ B2 2 HB
The quantity η B is also a regularization parameter. In analogy with the previous paragraph, the optimal solutions can be written as follows:

ψ Bj (•) = Nreg ℓ=1 γ ℓ,j κ B (W ℓ , •),
where κ B : IR V +3 × IR V +3 → IR is a reproducing kernel, and γ ℓ,j , ℓ ∈ {1, . . . , N reg } and j ∈ {1, 2}, are the unknown parameters to be determined. Let γ denote the N reg × 2 matrix whose (ℓ, j)-th entry is γ ℓ,j , and γ ℓ, * its ℓ-th row. Following the same line of reasoning as in the previous paragraph, one can collect the two estimations into a single one. Then, the unknown parameters γ ℓ,j are estimated at once, as follows:

γ = (K B + η B N reg I) -1 (X 0 Y 0 ), (15)
where

K B is the N reg × N reg matrix whose (v, w)-th entry is κ B (W v , W w ), for v, w ∈ {1, . . . , N reg }.
Finally, one is able to write ψ B as follows:

ψ B (•) = Nreg ℓ=1 γ ℓ, * κ B (W ℓ , •). (16)
After determining the first estimates of the parameters of the source s given by equation [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF], one can set the

(V + 3) × 1 vector W = (c (s) ⊤ Q Kx Kz ) ⊤ .
Using the model ψ B , a new estimate of the source's location is obtained as follows:

(x 0enh ŷ0enh ) = ψ B (W ), (17)
where (x 0enh ŷ0enh ) is the new enhanced location estimate of the source.

It is important to note that we are only using the first position estimate to compare its accuracy towards the enhanced one. When implementing the algorithm in a practical setup, one does not need to find the first position estimate, since it is not used in the enhancement process, which then reduces the complexity of the algorithm in terms of storage and time.

VI. SIMULATIONS AND RESULTS

In this section, we evaluate the performance of the proposed detection and estimation method for different scenarios. We consider a region of study of dimensions 5000 m × 5000 m. Sensors are deployed in a uniformly distributed manner in the area at a rate of one sensor each five meters, which is equivalent to a density of 0.05 sensor/m 2 . The region is divided into Z = 25 clusters with same size, thus having the same density of sensors. Fig. 4 shows the region's topology. The concentrations measured by the sensors are generated using the advection-diffusion model of [START_REF] Delmaire | Joint estimation decision methods for source localization and restoration in parametric convolution processes. application to accidental pollutant release[END_REF]. The maximum mass release Q max is taken equal to 1000 kg. The eddy diffusivity in the Z direction, K z is taken equal to 0.211 m 2 /s, while the eddy diffusivities K x and K y , in the X and Y directions respectively, are taken equal to 12 m 2 /s. As for the mean wind velocity U , it is taken equal to (1.8; 0; 0) m/s for illustrative purposes [START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF]. Finally, the sampling time T is taken equal to 1 s. The rest of this section is organized as follows. The training parameters are defined in the first subsection. In the second subsection, an evaluation of the method is proposed in the case of a single source with different noise values; then, in the third subsection, an evaluation in the case of multiple gas releases is provided. Finally, in the fourth subsection, a comparison of the proposed framework to the state-of-theart is provided.

A. Training parameters

First, we start by describing the training phase for the detection. As we already mentioned, we assume that the atmospheric conditions are the same in all the region of interest. Therefore, in order to reduce computations, we only define one detector; then, this detector is communicated to all cluster heads in order to be used in the clusters. The advantage is that only one training phase is needed now. The size of the training set for the detection N det is taken equal to 330, where 300 are normal data and 30 are outliers. We mean by outliers concentrations measured by sensors in the case of a gas release. In order to generate the outlier data, 30 different gas release scenarios are considered, where the source could be anywhere in the cluster, and the gas release mass Q is also randomly varying between 5 kg and Q max = 1000 kg. The Gaussian kernel is considered in our simulations. Consequently, the detection kernel κ : IR nz × IR nz → IR is given by the following:

κ(C i , C j ) = exp - C i -C j 2 2σ 2 ,
where i, j ∈ {1, ..., N det }, and σ is the kernel's bandwidth that plays a crucial role in defining the boundary around the training data.

Next, let us give a brief description of the training phase for the estimation step of the proposed framework. According to Section V, after finding the local maxima, only the concentrations around these maxima are used with the models ψ A and ψ B . By solving [START_REF] Mahfouz | Ridge regression and kalman filtering for target tracking in wireless sensor networks[END_REF], one can find the size of the zone encapsulating the useful concentrations. For Q max = 1000 kg and C thres = 10 -5 kg/m 3 , the solution is then ∆x = 26 m. Therefore, we only need to use a zone of 52 m × 52 m for the training. Let the number of considered scenarios N reg in the training phase be equal to 200, where x 0 ∈ [0; 52]m and y 0 ∈ [0; 52]m. Q is again randomly varying between 5 kg and 1000 kg. The Gaussian kernel is used here as well. The first kernel κ A : IR V ×IR V → IR is then given by the following:

κ A (c v , c w ) = exp - c v -c w 2 2σ 2 A ,
where v, w ∈ {1, ..., N reg }, and σ A is the kernel's bandwidth that controls, together with the regularization parameter η A , the degree of smoothness, noise tolerance, and generalization of the solution.

Due to the difference in the order of magnitude between the elements of the input W , we consider four different distances and bandwidths for the second kernel κ B : IR V +3 × IR V +3 → IR. Then, κ B is given by the following:

κ B (W v , W w ) = exp - cv -cw 2 2σ 2 B1 × exp - (Qv -Qw) 2 2σ 2 B2 × exp - (K xv -K xw ) 2 2σ 2 B3 × exp - (K z v -K z w ) 2 2σ 2

B4

, where v, w ∈ {1, ..., N reg }, and σ B1 , σ B2 , σ B3 , σ B4 are also the kernel's bandwidths, with σ B3 = σ B4 , since K x and K z have the same physical properties. The regularization parameters and the kernels' bandwidths are chosen in a way to minimize the error on the training set, using the 10-fold cross-validation technique [START_REF] Stone | Cross-validatory choice and assessment of statistical predictions[END_REF]. Simulation results show that taking σ B1 = σ B2 = σ B3 yields very close results to the ones obtained by taking different values for the bandwidths. Therefore, in order to save time and computations in the cross-validation phase, we take σ B1 = σ B2 = σ B3 .

B. Testing in the case of a single source

Now that we have defined all the parameters for the training of the detector and the two regression models, we can evaluate the performance of the proposed method in the case of a single source. The source is randomly positioned in the region of interest, with x 0 ∈ [0; 5000]m and y 0 ∈ [0; 5000]m, and Q ∈ [5; 1000] kg. All other parameters are set according to the beginning of this section, and as taken in [START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF]. The testing phase takes place during an interval of 10 s; then, the concentrations are measured each T = 1 s during this interval. The gas release can occur at any time in this interval. Note that random relative noises varying from 1% to 5% are added to the simulated concentration vectors, in order to account for transmission impairments or sensing errors. Using the proposed method, the gas release is detected in the affected clusters, and the source parameters are estimated based on the measured concentrations from the clusters. Table I shows the mean percentage of errors on the estimated source's parameters averaged over 50 Monte-Carlo simulations. We have a 100% detection in the affected zones, and one can see that the parameters of the source are accurately estimated. Also, notice the improvement in the source's location after the enhancement phase of Subsection V-B. Moreover, one can see from the small increase in the percentage of errors that the method is robust to noise.

C. Testing in the case of multiple sources

Now let us evaluate the accuracy of the method in the case of multiple gas releases, having different parameters and occurring at different times. To this end, consider four sources, whose parameters are given in Table II. The concentrations are generated using [START_REF] Delmaire | Joint estimation decision methods for source localization and restoration in parametric convolution processes. application to accidental pollutant release[END_REF], and a random relative noise of 5% is considered. The test phase is run here as well during an interval of 10 seconds. The subplots of Fig. 5 show the concentrations measured by the sensors around Source 2, whose parameters are given in Table II. Notice how the maximal concentration at each time t is moving along with the X-axis, and its value is decreasing with time because of the diffusion phenomenon. The estimated sources' parameters and time of detection of each source are given in Table III. One can see that the estimations are accurate, showing that the proposed method can also be used for the case of multiple sources.

D. Comparison to the state-of-the-art

In this subsection, we compare the percentage of errors to the ones obtained with the method in [START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF], for different values of the random relative noise varying from 1% to 5%. In [START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF], the source's parameters are estimated in an area of around 5500 m × 5500 m, using the measured concentrations and an inverse model obtained after solving a non-linear least squares estimation problem. The source is placed at a fixed location with x 0 = 5000 and y 0 = 100, and the release rate Q is taken equal to 1000 kg. Table IV shows the percentages of error for both methods, where PM denotes the proposed one and KM the one in [START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF]. Note that in this table, for the proposed method, only the enhanced estimates of the source's location are taken into account. Also, results over K z are not shown, since in [START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF], K z is assumed to be known. When comparing the errors, one can see that the proposed method greatly outperforms the KM method in terms of accuracy. Moreover, it is important to see here that the proposed method is much more robust to the increase of the noise percentage compared to the KM.

VII. CONCLUSION

This paper introduces a new clusterized method for the detection and estimation of the parameters of multiple gas sources in wireless sensor networks. Both phases of the method are defined within the framework of kernel methods, that proved to be very efficient in modeling nonlinear problems. Indeed, the evaluation of the proposed method on simulated data showed that the method yields accurate estimates, and the accuracy is maintained even with the increase of the noise level and in the case of multiple sources. Future work will handle further improvements of this method, such as to include cases where the eddy diffusivities or the wind's velocity and direction are not constant.

Fig. 1 :

 1 Fig. 1: Region topology-represents the cluster heads and • represents the deployed sensors.

Fig. 4 :

 4 Fig. 4: Region topology-represents the cluster heads.

Fig. 5 :

 5 Fig. 5: Concentration distribution in kg/m 3 at different times t around Source 2 (t 0 = 3 s).

TABLE I :

 I Percentage of errors on the estimated source's parameters for different random relative noises.

	Noise	Q	Kx	Kz	x0	ŷ0	x0enh ŷ0enh
	1%	0.69 0.006 0.006 0.02 0.03 0.01	0.01
	2%	0.75 0.007 0.007 0.03 0.04 0.01	0.02
	3%	0.89	0.01	0.01	0.03 0.04 0.02	0.02
	4%	0.96	0.02	0.02	0.04 0.04 0.02	0.02
	5%	1.15	0.02	0.02	0.04 0.04 0.02	0.02

TABLE II :

 II Parameters of the four sources.

		Source 1 Source 2 Source 3 Source 4
	t0 (s)	2	3	5	5
	Q (kg)	500	1000	20	200
	x0 (m)	100	100	2500	4500
	y0 (m)	100	4000	3000	1500
	Kx (m 2 /s)	12	12	12	12
	Kz (m 2 /s)	0.211	0.211	0.211	0.211

TABLE III :

 III Estimation of the parameters of the four sources using the proposed framework.

		Source 1 Source 2 Source 3 Source 4
	detection (s)	3	4	6	6
	Q (kg)	503.5	993	20.3	201.1
	x0enh (m)	100.0	100.2	2499.9	4500.3
	ŷ0enh (m)	100.1	3999.9	2999.2	1499.8
	Kx (m 2 /s)	12.01	11.99	11.98	11.99
	Kz (m 2 /s)	0.211	0.211	0.211	0.211

TABLE IV :

 IV Comparison of the percentage of errors obtained by the proposed method (PM) and the method in[START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF] (KM), for different random relative noises.

	Noise	Method	Q	Kx	x0	ŷ0
	1%	PM	0.32 0.01 0.01 0.02
		KM	2.7	0.7	0.8	0.4
	2%	PM	0.41 0.01 0.01 0.01
		KM	5.8	1.4	1.5	0.7
	3%	PM	0.52 0.01 0.01 0.02
		KM	8.0	2.2	2.4	1.1
	4%	PM	0.82 0.02 0.01 0.02
		KM	11.0	2.7	2.9	1.6
	5%	PM	1.05 0.02 0.01 0.03
		KM	12.3	3.7	4.0	1.9

This work is supported by the Champagne-Ardenne Region in France, grant WiDiD


Nisrine Ghadban

Paul Honeine

Farah Mourad-Chehade

Clovis Francis

Joumana Farah

In-network principal component analysis and diffusion strategies

Keywords: Principal component analysis, network, adaptive learning, distributed processing, dimensionality reduction

Principal component analysis is a very well-known statistical analysis technique. It determines principal axes that allow relevant data reduction. Thus, it is very useful in networks where a large number of measures is performed. However, such a technique needs high-complexity computation, which does not scales well in large networks. This paper proposes to overcome this issue, by describing several in-network strategies that estimate the principal axes with lower computational complexity. These strategies are either noncooperative or cooperative with information diffusion. The performance of the proposed strategies are illustrated on diverse applications, such as image processing and dimensionality reduction of time series in wireless sensor networks.

I. INTRODUCTION

Principal component analysis (PCA) is one of the most popular techniques for data analysis and processing, investigated in statistical analysis, data compression, and feature extraction [START_REF] Abdi | Principal component analysis[END_REF]. It consists in determining a subspace that retains the largest variance of the data. This subspace is spanned by the most relevant principal axes. Projecting the data onto the resulting subspace leads to the so-called principal components. By representing the data by these principal components, one reduces the dimensionality and extracts the data structure. The PCA technique is very useful in many applications involving multivariate analysis, e.g., data validation, fault detection [START_REF] Wise | Upset and sensor failure detection in multivariable processes[END_REF], and quality control [START_REF] Macgregor | Multivariate statistical methods for monitoring large data sets from chemical processes[END_REF]. It provides a powerful tool for data analysis and pattern recognition which is often used in signal and image processing [START_REF] Honeine | Online kernel principal component analysis: a reduced-order model[END_REF]. Moreover, in the context of sensor networks, the PCA has been investigated to extract features from noisy samples [START_REF] Chitradevi | Designing an efficient PCA based data model for wireless sensor networks[END_REF], compress and denoise time series measurements [START_REF] Chen | Algorithm of data compression based on multiple principal component analysis over the wsn[END_REF], [START_REF] Rooshenas | Reducing the data transmission in wireless sensor networks using the principal component analysis[END_REF], as well as for detection of intrusion [START_REF] Ahmadi Livani | A PCA-based distributed approach for intrusion detection in wireless sensor networks[END_REF] or anomaly [START_REF] Huang | In-network pca and anomaly detection[END_REF].

The conventional PCA requires the eigen-decomposition of the sample covariance matrix, where the eigenvectors associated to the largest eigenvalues are the most relevant principal axes. Therefore, it is necessary to send all data to a fusion center (FC), where the covariance matrix is estimated and eigen-decomposed. Such configuration is not scalable and inappropriate for network applications. Moreover, the computational complexity of such eigendecomposition is cubic with the size of the dataset. Several attempts have been made to alleviate this problem, such as in [START_REF] Bunch | Updating the singular value decomposition[END_REF], [START_REF] Hall | Merging and splitting eigenspace models[END_REF] where algorithms are investigated to compute the eigenspace, but still with high computational cost. More recently, in [START_REF] Kargupta | Collective principal component analysis from distributed, heterogeneous data[END_REF], a technique was proposed in which the FC only receives the principal components, instead of the whole time series. In [START_REF] Le Borgne | Distributed principal component analysis for wireless sensor networks[END_REF], the most relevant principal axis is estimated by the power iteration method. This method requires the computation of the sample covariance matrix, making it inappropriate for in-network processing.

In this paper, we propose to estimate the most relevant principal axes without computing the sample covariance matrix. To this end, we investigate the Oja's neural-based rule [START_REF] Oja | Simplified neuron model as a principal component analyzer[END_REF], which has been recently revisited in [START_REF] Honeine | Online kernel principal component analysis: a reduced-order model[END_REF] for nonlinear PCA with kernel-based machines. We extend our study in [START_REF] Ghadban | Strategies for principal component analysis in wireless sensor networks[END_REF] to incorporate several cost functions, including the information theoretical criterion derived in [START_REF] Plumbley | Lyapunov functions for convergence of principal component algorithms[END_REF], Rayleigh quotient, LUO and OJAn cost functions derived in [START_REF] Cirrincione | The mca exin neuron for the minor component analysis[END_REF], [START_REF] Luo | A minor component analysis algorithm[END_REF]. Within a network setting, several in-network strategies are proposed, including noncooperative fusion center

Centralized network

Noncooperative network Cooperative (diffusion) network Fig. 1. Illustration of information processing in a network: While the centralized network requires a fusion center, the noncooperative and cooperative networks provide in-network processing. The noncooperative network process the information by using a given routing path.

The cooperation network operates information diffusion, using either a combine-then-adapt or an adapt-then-combine strategy.

and cooperative strategies. In the latter, sensors cooperate with each other by information diffusion in order to estimate the principal axes. We derive two diffusion strategies, the so-called combine-then-adapt and adapt-thencombine strategies, which have been recently investigated in linear adaptive filtering literature by A. Sayed et al.

in [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks: an examination of distributed strategies and network behavior[END_REF]. To the best of our knowledge, our study is the first work that investigates these adaptation strategies for unsupervised learning, as given here with the PCA. To provide a comprehensive presentation throughout this paper, we first study in detail the estimation of a single principal axis for a given cost function, before generalizing it to different cost functions. We then extend the proposed approach to the estimation of multiple principal axes, by investigating orthogonalization processes, such as the well-known Gram-Schmidt orthogonalization and the deflation scheme.

The rest of this paper is organized as follows. Section II briefly introduces the network topology. We describe the proposed strategies for estimating the first principal axis based on Oja's rule in Section III. In Section IV, we generalize the proposed approach to other cost functions, beyond Oja's rule. Section V extends the strategies for multiple axes extraction. In Section VI, we present algorithms well-known in the literature. Section VII provides experimentation results and discussions, whereas Section VIII concludes the paper.

II. NETWORK SETTING

Consider a network of N connected nodes (also called agents in the literature). In such a network, any two nodes are necessarily linked: either directly if they are neighbors, or through other intermediate nodes. On the other hand, there are mainly two types of networks: centralized and decentralized networks. In the former, the nodes are connected to a FC, and in the latter the nodes are connected in a noncooperative way by a given routing path, or in a cooperative way, where each node communicates with its neighbors. We denote by V k the set of indices of the neighboring nodes to node k, with k = 1, ..., N , i.e., the nodes that are directly connected to it. We consider that the node k is adjacent to itself, that is to say k ∈ V k . The network topologies studied in this paper are illustrated in Fig. 1.

We are interested in the case where the nodes estimate the same principal axis w based on some measured information corresponding to a common phenomenon. Let x k be the p-by-1 vector collected by the node k, for k = 1, . . . , N , p being the measurements' dimension, and let X ⊂ R p be the space of these collected data (assumed to be zero-mean), with the conventional inner product x ⊤ k x l for any x k , x l ∈ X. We denote by the scalar value y w,k = w ⊤ x k the inner product associated with the orthogonal projection of any x k ∈ X onto the vector w ∈ X and by y w the scalar random variable taking the values y w,k , with k = 1, ..., N . The ultimate goal would be to compute w, optimaly for some given cost function and according to the network topologies, to keep afterwards only y w,k and w, for k = 1, ..., N , allowing us to represent efficiently the data x k when needed.

III. IN-NETWORK PCA

In this section, we describe the different strategies depending on the network topology. We first consider the centralized network and expose the conventional strategy for PCA. We then derive our new strategies applied to the decentralized networks for extracting the first principal axis. We extend these techniques to the case of multiple principal axes in Section V.

A. Centralized strategy

Here, the nodes are assumed to be connected to a FC. The latter receives the collected data without any local processing from the nodes. In order to extract the first principal axis, the FC maximizes the variance of the projected data (assumed to be zero-mean), namely max w E(y 2 w ), where E(•) is the expectation over the density of the input data. By taking the empirical estimation of the variance of the projected data with respect to the available samples, x 1 , x 2 , . . . , x N , we get max w w ⊤ Cw, where

C = 1 N N k=1 x k x ⊤
k is the covariance of the measured data. The problem takes the form Cw = λw. This is the well-known eigen-decomposition problem where w is the eigenvector associated with the eigenvalue λ of C. In fact, the latter is the corresponding variance of the projected data, since

w ⊤ Cw = w ⊤ λw = λw ⊤ w = λ.
Therefore, we must consider the eigenvector associated with the largest eigenvalue to maximize this variance. Note that the k-th principal component of x is w ⊤ k x and var[w ⊤ k x] = λ k , where the k-th principal axis w k is the eigenvector corresponding to the k-th largest eigenvalue λ k .

Having the data, x 1 , . . . , x N , the FC solves the eigen-decomposition problem. The eigenvector associated with the largest eigenvalue corresponds to the principal axis, denoted w * in the following. The computational complexity of the eigen-decomposition problem is O(p 3 ). There is also the communication complexity, which is O(N p) over a distance O(1). In the following, we propose strategies that avoid the computation of the sample covariance matrix and its eigen-decomposition, thus allowing to reduce significantly the computational complexity.

B. Noncooperative strategy

We consider an in-network scheme where the first principal axis is learned adaptively instead of directly maximizing the overall projected variance E(y 2 w ) and explicitly investigating the covariance matrix. To this end, we consider an "instantaneous" estimation of the principal axis w at each node. According to a given routing process, each node k receives an estimate w t-1 from another node, and adjusts it using its own data x k by maximizing y 2 w,k . In other words, the node k minimizes the "instantaneous" quadratic reconstruction error, namely

J k (w) = 1 4 x k -y w,k w 2 . ( 1
)
The gradient of this cost function with respect to w is

∇ w J k (w) = y w,k (y w,k w -x k ), (2)
which leads to the gradient descent formulation

w t = w t-1 + η t (x k y wt-1,k -y 2 wt-1,k w t-1 ), (3)
where η t is the learning rate. It turns out that this update rule is essentially Oja's rule [START_REF] Oja | Simplified neuron model as a principal component analyzer[END_REF], which is a single-neuron special case of the generalized Hebbian learning [START_REF] Dony | Neural network approaches to image compression[END_REF], applied here within the in-network setting described above. This learning rule converges to an equilibrium state, which is the first principal axis w * . In fact, when w t converges to some state w, we have from (3) that x k y w,k = y 2 w,k w, namely x k x ⊤ k w = y 2 w,k w. Averaging the latter expression over the whole data, we get the well-known eigen-decomposition problem of the covariance matrix Cw = E(y 2 w ) w, where the eigenvalue is the squared output y w,k to be maximized. Therefore, the update rule (3) converges to the eigenvector associated with the largest eigenvalue of the covariance matrix.

C. Cooperative strategies

We consider now a cooperative strategy where each node k has access to the information of its neighborhood V k , which is the subset of nodes directly connected to node k, including itself. The optimization problem can be presented as a minimization of the cost function N k=1 J k (w), where J k (w) is the cost function at node k, for instance [START_REF] Abdi | Principal component analysis[END_REF]. Since each node k communicates only with its neighboring nodes, we introduce the nonnegative coefficients c kl that relate node k to its neighbors. They satisfy the following conditions:

c kl ≥ 0, l∈Vk c kl = 1,
and

c kl = 0 if l / ∈ V k . (4)
It is worth noting that if l / ∈ V k , k / ∈ V l and thus c lk is also null. Using these coefficients, we define, for each node l, a local cost that consists of a weighted combination of the individual costs of the neighbors of node l (including l itself):

J loc l (w) = k∈Vl c kl J k (w) = N k=1 c kl J k (w). (5)
It turns out that the cumulative sum of these local cost functions is equal to the global cost function, since we have

N k=1 J k (w) = N k=1 N l=1 c kl J k (w) = N l=1 N k=1 c kl J k (w) = N l=1 J loc l (w).
For node k, the global cost function can be decomposed as

N l=1 J l (w) = J loc k (w) + N l=1 l =k J loc l (w). (6)
While the first term in the right-hand-side is known at the node under scrutiny, the second one should be estimated using information from the neighbors of node k; thus, the above summation is restricted to its neighborhood. Moreover, we relax it by constraining the norm between the estimated vector w and the optimal (global) principal axis w * within the neighborhood of node k, with the following regularization:

N l=1 l =k J loc l (w) ≈ l∈Vk\{k} b lk w -w * 2 ,
where the parameters b lk control the tradeoff between the accuracy and the smoothness of the solution. Such regularization is also motivated by investigating the second-order Taylor expansion of J loc l (w), as described in [START_REF] Chen | Diffusion adaptation strategies for distributed optimization and learning over networks[END_REF]. Note that we use w * knowing that we do not have access to its value. We shall show in the following how to overcome this problem.

By injecting this approximation into (6), and using ( 5), the minimization of the global cost function at k is equivalent to minimizing

J glob k (w) = J loc k (w) + l∈Vk\{k} b lk w -w * 2 = l∈Vk c lk J l (w) + l∈Vk\{k} b lk w -w * 2 , (7)
where the last equality follows from [START_REF] Chitradevi | Designing an efficient PCA based data model for wireless sensor networks[END_REF]. In order to minimize the above cost function, the node k applies the gradient descent to J glob k (w) with:

w k,t = w k,t-1 -η k,t ∇ w J glob k (w k,t-1 ). (8)
Here, η k,t is the learning rate for node k at iteration t. Replacing J glob k (w) by its expression in [START_REF] Rooshenas | Reducing the data transmission in wireless sensor networks using the principal component analysis[END_REF], we get:

w k,t = w k,t-1 -η k,t l∈Vk c lk ∇ w J l (w t-1 ) + η k,t l∈Vk\{k} b lk (w * -w k,t-1 ).
In this expression, ∇ w J l (w t-1 ) is the gradient of the PCA-based cost function, presented earlier in [START_REF] Wise | Upset and sensor failure detection in multivariable processes[END_REF]. In this case, we get

w k,t = w k,t-1 + η k,t l∈Vk c lk x l y wt-1,l -y 2 wt-1,l w l,t-1 + η k,t l∈Vk\{k} b lk (w * -w k,t-1 ). (9)
Without loss of generality, we shall consider the update rule (9) in the following. This update rule from w k,t-1 to w k,t involves adding two correction terms: the first one operates in the maximum variance direction and the second one associates neighborhood regularization. By decomposing this calculation into two successive steps, we get two possible strategies, depending on the order of adding the correction terms, and that differ essentially in the approximation of the unknown w * in (9), as shown next.

Adapt-then-combine (ATC) strategy

We express the update rule (9) as follows:

φ k,t = w k,t-1 + η k,t l∈Vk c lk x l y wt-1,l -y 2 wt-1,l w l,t-1 , w k,t = φ k,t + η k,t l∈Vk\{k} b lk (w * -w k,t-1 ).
The first step uses local gradient vectors from the neighborhood of the node k in order to update w k,t-1 to the intermediate estimate φ k,t . The second step updates φ k,t to w k,t . This second step is not realizable since the nodes do not know w * . However, each node l has its own approximation of w * , which is its local intermediate estimate φ l,t . Therefore, we replace w * by φ l,t . On the other hand, the intermediate value φ k,t at node k is obtained by incorporating information from the neighbors (as given by the first step). Thus, it is generally a better estimate for w * than w k,t-1 . Therefore, we further replace w k,t-1 by φ k,t in the second step. Hence, the second step is rewritten as:

w k,t = φ k,t + η k,t l∈Vk\{k} b lk (φ l,t -φ k,t ) = 1 -η k,t l∈Vk\{k} b lk φ k,t + η k,t l∈Vk\{k} b lk φ l,t ,
or, equivalently: w k,t = l∈Vk a kl φ l,t where we have used the following nonnegative weighting coefficients:

a kl =    1 -η k,t i∈Vk\{k} b ik , if l = k; η k,t b lk , if l ∈ V k \{k} ; 0, otherwise.
Finally, the update rule for the adapt-then-combine strategy is:

φ k,t = w k,t-1 + η k,t l∈Vk c lk x l y wt-1,l -y 2 wt-1,l w l,t-1 , w k,t = l∈Vk a kl φ l,t .
In the special case where no information exchange is performed between nodes, except for the aggregation step, i.e., c kk = 1 and c kl = 0 for all l = k, we get:

φ k,t = w k,t-1 + η k,t x k y wt-1,k -y 2 wt-1,k w k,t-1 w k,t = l∈Vk a kl φ l,t .
Note that such particular configuration would allow minimizing the amount of information exchanged between nodes, by excluding the transmission of the estimates y and the measurements x within the network, and restricting the exchanges to w (or equivalently φ).

Combine-then-adapt (CTA) strategy

In this strategy, we express the update rule (9) as follows:

φ k,t = w k,t-1 + η k,t l∈Vk\{k} b lk (w * -w k,t-1 ) w k,t = φ k,t + η k,t l∈Vk c lk x l y w t-1 ,l -y 2 wt-1,l w l,t-1 .
For the same reasons shown in the adapt-then-combine strategy, we replace in the first step w * by w l,t-1 and, in the second step, we replace w l,t-1 by φ l,t . Introducing the same coefficients a kl , we obtain the update rule for the combine-then-adapt strategy:

φ k,t = l∈Vk a kl w l,t-1 , w k,t = φ k,t + η k,t l∈Vk c lk x l y wt-1,l -y 2 wt-1,l φ l,t .
In the case of no information exchange except for the aggregation step, we obtain:

φ k,t = l∈Vk a kl w l,t-1 w k,t = φ k,t + η k,t x k y w t-1 ,k -y 2 wt-1,k φ k,t .

D. Connections to the consensus strategies

The consensus-type implementation [START_REF] Degroot | Reaching a Consensus[END_REF] is a class of distributed strategies that reveals the coordination of networks of autonomous agents. In order to estimate a common parameter w * , the consensus allows each node k to align its decision w k,t with the decisions of its neighbors:

w k,t = l∈Vk a kl w k,t-1 ,
where a kl are nonnegative weights that node k gives to the estimates received from its neighbors. In order to distribute the computations among the nodes, the authors in [START_REF] Tu | Diffusion strategies outperform consensus strategies for distributed estimation over adaptive networks[END_REF], [START_REF] Nedic | Cooperative distributed multi-agent optimization[END_REF], [START_REF] Yuan | On the convergence of decentralized gradient descent[END_REF] propose to combine the consensus and the gradient descent approaches. In such algorithm, each node k applies the negative gradient step to the average of its neighbors' estimates. Therefore, the node k aligns its decision w k,t with the decisions of its neighbors and minimizes its own objective J k (w) using:

w k,t = l∈Vk a kl w k,t-1 -η k,t ∇ w J k (w k,t-1 ).
In this context, we revisit the estimation of the principal axis with the proposed cost function given in (1). This strategy takes the following form:

w k,t = l∈Vk a kl w k,t-1 + η k,t x k y wt-1,k -y 2 wt-1,k w k,t-1 .
The consensus strategy, the adapt-then-combine, and the combine-then-adapt diffusion strategies have the same computational complexity. However, the diffusion strategies outperform the consensus implementation. Indeed, the adapt-then-combine and the combine-then-adapt diffusion strategies can be expressed respectively by:

w k,t = l∈Vk a kl w k,t-1 + η k,t x k y w,k -y 2 w,k w k,t-1 , w k,t = l∈Vk a kl w k,t-1 + η k,t x k y w,k -y 2 w,k l∈Vk a kl w k,t-1 .
Considering for instance the combine-then-adapt strategy, we note that the combination coefficients a kl appear in the correction term, while the consensus uses only w k,t-1 to correct the error (i.e., without involving the neighborhood), which is less effective, as shown in the experimental results.

IV. OTHER COST FUNCTIONS

In this part, we extend the study of Section III to include different cost functions. These cost functions are chosen in the light of some known rules, as described next.

A. Information theory (IT)

Drawing inspiration from the information theoretical framework studied in [START_REF] Plumbley | Lyapunov functions for convergence of principal component algorithms[END_REF], [START_REF] Miao | Fast subspace tracking and neural network learning by a novel information criterion[END_REF], we consider the maximization of the entropy, defined by ln(y 2 w,k ), under a constraint on the norm of w. To this end, we minimize the following cost function associated with node k:

J k (w) = 1 2 w ⊤ w -1 2 ln(y 2 w,k ), (10)
whose gradient with respect to w is

∇ w J k (w) = w - x k y w,k . (11)
By applying the gradient descent technique to this cost function, we obtain the following noncooperative update rule:

w t = w t-1 + η t x k y wt-1,k -w t-1 .
Cooperative strategies are also investigated by information diffusion. The adapt-then-combine update rule becomes in this case:

φ k,t = w k,t-1 + η k,t l∈Vk c lk x l y w t-1 ,l -w l,t-1 , w k,t = l∈Vk a kl φ l,t .
The combine-then-adapt update rule becomes:

φ k,t = l∈Vk a kl w l,t-1 , w k,t = φ k,t + η k,t l∈Vk c lk x l y wt-1,l -w l,t-1 .
While the cost function (1) studied in Section III and the cost function [START_REF] Kargupta | Collective principal component analysis from distributed, heterogeneous data[END_REF] inspired from the information theory are very different, the corresponding update rules are intimately connected. Indeed, the latter can be derived from the former by using the normalized learning rate η t /y 2 wt-1,k .

B. Rayleigh quotient (RQ)

The eigen-decomposition problem is highly connected to the Rayleigh quotient. Inspired from the Rayleigh quotient framework studied in [START_REF] Cirrincione | The mca exin neuron for the minor component analysis[END_REF], [START_REF] Luo | A minor component analysis algorithm[END_REF], we minimize the following cost function associated with node k:

J k (w) = - y 2 w,k 2 w ⊤ w , (12)
whose gradient with respect to w is

∇ w J k (w) = 1 w ⊤ w w y 2 w,k w ⊤ w -x k y w,k . (13)
The gradient descent technique, applied to this cost function, leads to the following noncooperative update rule:

w t = w t-1 + η t 1 w ⊤ t-1 w t-1 x k y wt-1,k -w t-1 y 2 wt-1,k w ⊤ t-1 w t-1
.

By including a diffusion process, we propose two cooperative strategies. In the case of the adapt-then-combine update rule, we get:

φ k,t = w k,t-1 + η k,t l∈Vk c lk 1 w ⊤ l,t-1 w l,t-1 x l y wt-1,l -w l,t-1 y 2 wt-1,l w ⊤ l,t-1 w l,t-1 , w k,t = l∈Vk a kl φ l,t .
For the combine-then-adapt update rule, we have:

φ k,t = l∈Vk a kl w l,t-1 , w k,t = φ k,t + η k,t l∈Vk c lk 1 w ⊤ l,t-1 w l,t-1 x l y wt-1,l -w l,t-1 y 2 wt-1,l w ⊤ l,t-1 w l,t-1
.

C. OJAn and LUO algorithms

The OJAn and LUO algorithms, presented by Luo et al. in [START_REF] Luo | A minor component analysis algorithm[END_REF], are inter-connected to the Rayleigh quotient [START_REF] Cirrincione | The mca exin neuron for the minor component analysis[END_REF], [START_REF] Taleb | Against the convergence of the minor component analysis neurons[END_REF]. The OJAn algorithm is a normalized version of Oja's rule. It can be obtained from the Rayleigh quotient by using a learning rate η t (w ⊤ t-1 w t-1 ). By revisiting the gradient descent formulation (3) in the light of this learning rate, we get the following noncooperative update rule:

w t = w t-1 + η t x k y wt-1,k -w t-1 y 2 wt-1,k w ⊤ t-1 w t-1
.

In this case, the adapt-then-combine update rule becomes:

φ k,t = w k,t-1 + η k,t l∈Vk c lk x l y wt-1,l -w l,t-1 y 2 wt-1,l w ⊤ l,t-1 w l,t-1 , w k,t = l∈Vk a kl φ l,t .
The combine-then-adapt update rule becomes:

φ k,t = l∈Vk a kl w l,t-1 , w k,t = φ k,t + η k,t l∈Vk c lk x l y w t-1 ,l -w l,t-1 y 2 wt-1,l w ⊤ l,t-1 w l,t-1
.

The LUO algorithm can be derived from the Rayleigh quotient by using the learning rate η t (w ⊤ t-1 w t-1 ) 2 . This leads to following noncooperative update rule:

w t = w t-1 + η t (w ⊤ t-1 w t-1 ) x k y wt-1,k -w t-1 y 2 wt-1,k w ⊤ t-1 w t-1
, as well as the two cooperative strategies, the adapt-then-combine

φ k,t = w k,t-1 + η k,t l∈Vk c lk (w ⊤ l,t-1 w l,t-1 ) x l y wt-1,l -w l,t-1 y 2 wt-1,l w ⊤ l,t-1 w l,t-1 , w k,t = l∈Vk a kl φ l,t .
and the combine-then-adapt update rule

φ k,t = l∈Vk a kl w l,t-1 , w k,t = φ k,t + η k,t l∈Vk c lk (w ⊤ l,t-1 w l,t-1 ) x l y wt-1,l -w l,t-1 y 2 wt-1,l w ⊤ l,t-1 w l,t-1
.

V. EXTRACTING MULTIPLE PRINCIPAL AXES

This section extends the previous derivations in Section III to multiple principal axes, generalization to other cost functions as given in Section IV is straightforward. To this end, we denote by

W t = [w 1,t w 2,t • • • w r,t ] ⊤
the r-by-p matrix of the first r principal axes estimated at iteration t, sorted in descending order of their eigenvalues. Let

y W,k = [y w1,k y w2,k • • • y wr,k ] ⊤ ,
where y w j ,k = w ⊤ j,t x k . Next, we combine the update rules given in Sections III-B and III-C with several orthogonalization processes, including the well-known Gram-Schmidt orthogonalization and the deflation [START_REF] Möller | First-order approximation of gram-schmidt orthonormalization beats deflation in coupled pca learning rules[END_REF]. Note that other orthogonalization processes exist in the literature, such as the symmetric orthogonalization [START_REF] Srivastava | A unified view of the orthogonalization methods[END_REF]. However, the latter is inappropriate in network applications due to its high computational complexity.

A. Gram-Schmidt orthogonalization process

We use a generalized Hebbian approach as proposed by Sanger for the linear principal component analysis [START_REF] Sanger | Optimal unsupervised learning in a single-layer linear feedforward neural network[END_REF], [START_REF]Two iterative algorithms for computing the singular value decomposition from input/output samples[END_REF]. In the noncooperative strategy, we write the update rule of the j-th principal axis as follows:

w j,t = w j,t-1 + η t x k y wj,k -y wj,k j l=1 y wj,l w l,t-1 . (14)
It is clear how the estimation of the j-th component w j,t at iteration t involves the previous estimates (at iteration t -1) of lower-order components, i.e., w l,t-1 , for l = 1, . . . , j. By collecting all these estimates in a single matrix W t , we obtain the following update rule:

W t = W t-1 + η t y W t-1 ,k x ⊤ k -LT(y W t-1 ,k y ⊤ Wt-1,k ) W t-1 , (15)
where LT(•) makes its argument lower triangular by setting to zero the entries above its diagonal. Note that the learning rate does not need to be the same for all the principal axes.

In the cooperative strategies, we denote by

Φ k,t = [φ 1,k,t φ 2,k,t • • • φ r,k,t
] ⊤ the r-by-p matrix of the r intermediate estimates. In matrix form, the update rule associated with the adapt-then-combine strategy is as follows

Φ k,t = W k,t-1 + η k,t y Wt-1,k x ⊤ k -LT(y Wt-1,k y ⊤ Wt-1,k ) W k,t-1 , W k,t = l∈Vk a kl Φ l,t .
The update rule for the combine-then-adapt strategy becomes

Φ k,t = l∈Vk a kl W l,t-1 , W k,kt = Φ k,t + η k,t y Wt-1,k x ⊤ k -LT(y Wt-1,k y ⊤ Wt-1,k ) Φ k,t .

B. Deflation scheme

Like the Gram-Schmidt orthogonalization, the deflation scheme updates all estimates at each iteration. However, it proceeds in two consecutive steps, within each iteration. First, intermediate estimates w j,⋆ , for j = 1, . . . , r are independently calculated using:

w j,⋆ = w j,t-1 + η t (x k y wj,k -y 2 wj,k w j,t-1 ). (16)
Let

W ⋆ = [w 1,⋆ w 2,⋆ • • • w r,⋆
] be the matrix of the intermediate estimates. We can write:

W ⋆ = W t-1 + η t x k y ⊤ W,k -W t-1 diag(y W,k y ⊤ W,k ) ,
where diag(•) makes its argument diagonal by setting all non-diagonal entries to zero. Algorithm Time constant for w 1 Time constant for

w i Oja 1/λ 1 1/(λ 1 -λ i ) IT 1 1 RQ w k,0 2 /λ 1 w k,0 2 /(λ 1 -λ i ) OJAn 1/λ 1 1/(λ 1 -λ i ) LUO w k,0 -2 /λ 1 w k,0 -2 /(λ 1 -λ i )
Then, each estimate is updated by the node k using the intermediate estimates from [START_REF] Plumbley | Lyapunov functions for convergence of principal component algorithms[END_REF] as follows:

w j,t = w j,⋆ - j-1 l=1 (w ⊤ l,t-1 w j,⋆ ) w l,t-1 .
By regrouping all these vectors w j,t in a matrix W t , this expression for noncooperative strategy is written in matrix form as:

W t = W ⋆ -W t-1 W ⊤ t-1 W ⋆ -LT(W ⊤ t-1 W ⋆ ) .
A normalization of w j,t is required by the orthogonalization procedure, by operating the projection wj,t wj,t . Note that this scheme could be applied in a sequential way: The node k can extract the principal axes one after the other, by operating on a simple vector update. For the j-th principal axis, for j = 1, . . . , r, the intermediate j-th estimate defined in ( 16) is used. Then, the node under scrutiny updates the j-th principal axis with the following update rule:

w j,t = w ⋆ - j-1 l=1 (w ⊤ l,∞ w j,⋆ ) w l,∞ , (17)
where w l,∞ is the extracted l-th principal axis obtained after convergence of w l,t , for l = 1, . . . , j -1. However, such process is ill-suited for real-time network applications, because the extraction of the (j + 1)-th principal axis cannot start before the convergence of w j,t .

For the cooperative strategies, the deflation scheme can be easily incorporated in the adapt-then-combine strategy, with

   Φ k,⋆ = W k,t-1 + η k,t x k y ⊤ Wt-1,k -W k,t-1 diag(y W t-1 ,k y ⊤ Wt-1,k ) Φ k,t = Φ ⋆ -Φ k,t-1 Φ ⊤ k,t-1 Φ k,⋆ -LT(Φ ⊤ k,t-1 Φ k,⋆ ) W k,t = l∈Vk a kl Φ l,t .
as well as the combine-then-adapt strategy, with

Φ k,t-1 = l∈Vk a kl W l,t-1    W k,⋆ = Φ k,t-1 + η k,t x k y ⊤ Wt-1,k -Φ k,t-1 diag(y W t-1 ,k y ⊤ Wt-1,k ) W k,t = W ⋆ -W k,t-1 W ⊤ k,t-1 W k,⋆ -LT(W ⊤ k,t-1 W k,⋆ ) .
criterion for comparison is the level of the floor. The IT algorithm has the largest average angle at the convergence.

All the other algorithms have comparable levels of the floor. Note that the consensus with RQ, OJAn, and LUO has a smaller average angle at the convergence than the other strategies, and it is slower in convergence than the other strategies with Oja, OJAn and LUO. Altogether, we can say that the best combinations that allow for a small angle floor, together with an acceptable rate of convergence, are the Oja and OJAn rules with the adapt-then-combine diffusion strategy.

B. Image processing application

In this section, we study an image processing application. For this purpose, we consider a dataset of the handwritten digit "1" given by images of 28-by-28 pixels, treated as 784-dimensional vectors. In order to compare with the "ground truth" solution obtained from the eigen-decomposition of the sample covariance matrix, the number of images is restricted to 90. This time, the connectivity between nodes is tested for s = 2240. For the stepsize parameters, we take η 1 = 5 • 10 -8 for the first principal axis and a different one, η 2 = 1 • 10 -7 for the second principal axis, in order to insure a better convergence. Note that η 2 ≥ η 1 because the estimation of the second principal axis w 2,t requires the estimation of the first principal axis w 1,t-1 . The same stepsize values are taken for cooperative and noncooperative strategies. As for the IT algorithm, we take η 1 = η 2 = 0.05, for the Rayleigh quotient algorithm η 1 = η 2 = 5 • 10 -6 , for the OJAn algorithm η 1 = η 2 = 5 • 10 -8 , and for the LUO algorithm η 1 = 5 • 10 -11 and η 2 = 2.5 • 10 -10 , in cooperative and noncooperative strategies.

Fig. 4 shows the convergence of the different orthogonalization processes using Oja's rule. It shows that the Gram-Schmidt orthogonalization and the deflation schemes have a comparable efficiency, especially with the cooperative strategies. The Gram-Schmidt orthogonalization is slightly better than the deflation scheme. In what follows, we use the Gram-Schmidt orthogonalization for multiple axes.

Fig. 5 to Fig. 7 show the convergence of each strategy for Oja's rule and for the IT, the Rayleigh quotient, the OJAn, and the LUO algorithms respectively, with the angles averaged over the N nodes. Again, the noncooperative strategy is significantly outperformed by all diffusion strategies. On the other hand, considering the rate of convergence, the figures show that the OJAn, the Oja, and the IT algorithms have the fastest convergence, followed consecutively by the RQ and the LUO algorithms. In terms of the floor level, the IT algorithm has the largest average angle at the convergence. All the other algorithms have a comparable level of the floor. As in the WSN application given in Section VII-A, the combination of the Oja or the OJAn algorithms with the adapt-then-combine diffusion strategy seems to present the best compromise between the angle floor and the rate of convergence.

VIII. CONCLUSION

In this paper, we proposed several strategies for estimating the principal axes for PCA in networks. Noncooperative and cooperative strategies with information diffusion were studied for this purpose. Experiments were conducted by taking into account the constraints imposed in WSNs and image processing applications. The results showed the relevance of these strategies. The convergence of all proposed algorithms was better than the non-cooperative strategy, in the case of WSNs. The main reason for this result is the ratio of the first and second eigenvalues which is greater in the case of wireless sensor networks. Furthermore, the level of efficiency of certain cost functions, like with the RQ and the LUO, highly depends on the choice of the initialization of the principal axes, which is not the case for the Oja and the OJAn rules. A similar behavior was observed for the image processing application. As a future work, we will also include the use of the spatial information within the study.

IX. ACKNOWLEDGEMENT

This work is supported by the Région Champagne-Ardenne (grant "WiDiD") and the Lebanese University.

Fig. 2 .

 2 Fig. 2. Convergence analysis for Oja's rule, in the WSN settings.

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Convergence analysis in the WSN settings of the strategies based on information theory (top left), on Rayleigh quotient (top right), on OJAn (bottom left) and on LUO (bottom right) algorithms.

Fig. 6 .

 6 Fig. 6. Convergence analysis of the strategies based on the information theory (top figure) and the Rayleigh quotient (bottom figure), for the image dataset.

Fig. 7 .

 7 Fig. 7. Convergence analysis of the strategies based on the OJAn algorithm (top figure) and the LUO algorithm (bottom figure), for the image dataset.

TABLE I EXPRESSIONS

 I OF THE TIME CONSTANTS FOR THE FIRST AND THE i-TH PRINCIPAL AXES, FOR EACH OF THE ALGORITHMS STUDIED IN THIS PAPER.

Remarks on the convergence

The previous study takes into account the cost function [START_REF] Abdi | Principal component analysis[END_REF] and can be easily generalized to other algorithms by considering the cost function. In the light of the work in [START_REF] Chatterjee | Adaptive algorithms for first principal eigenvector computation[END_REF], it is worth noting some observations on the convergence of these algorithms. We refer the reader to [START_REF] Chatterjee | Adaptive algorithms for first principal eigenvector computation[END_REF] and references therein for more details.

The convergence of the IT algorithm is independent of the largest eigenvalue λ 1 . The convergence of all other algorithms increases with λ 1 and λ 1 /λ 2 , when extracting the first and second principal axes. Moreover, the convergence of Oja's based rule and the OJAn algorithm cannot be improved by increasing w k,0 , while the convergence of the LUO algorithm decreases for smaller w k,0 and the convergence of the RQ algorithm decreases for larger w k,0 . TABLE I shows the the rate of convergence in terms of the expressions of the time constants for the first and the i-th principal axes, for all the algorithms studied in this paper.

The algorithms associated to Oja's rule and to the IT converge to unit norm vectors, i.e., w k,∞ → 1. On the other hand, the RQ, OJAn, and LUO algorithms require a normalization w k,0 = 1 for the principal eigenvector w k,∞ to converge to a unit-norm solution.

VI. COMPARISON WITH OTHER ALGORITHMS

In this section, we present two algorithms well-known in the literature, collective-PCA and PCA-based distributed approach. As opposed to these algorithms where one or many eigen-decomposition problems are solved, the innetwork noncooperative and cooperative strategies proposed in this paper neither estimate the covariance matrix nor operate an eigen-decomposition. Experimental results are given in the following section.

A. Collective-PCA

The collective-PCA (CPCA) introduced in [START_REF] Kargupta | Collective principal component analysis from distributed, heterogeneous data[END_REF] essentially operates in two stages: pre-processing the data in cluster-heads (i.e., local fusion centers) before sending them to a FC. Let X be the N -by-p matrix of the entire data, i.e., X =

where X i is a N -by-p i submatrix of X, and d i=1 p i = p. A local centralized PCA is performed on the data of each X i . Let W i be the p i -by-r i matrix of the r i first principal axes, and W tot be a diagonal matix with matrices W 1 , . . . , W d on the diagonal. The principal component corresponding to X

The FC performs a PCA on Y and extract its principal axes v i . Finally, the principal axes of X are w i = W tot v i . This centralized algorithm reduces the complexity of the communication with the FC. However, it still has an important computational complexity compared to our algorithms, because it requires to perform centralized PCA on different stages.

B. PCA-based distributed approach

The PCA-based distributed approach (PCADID) proposed in [START_REF] Livani | A PCA-based distributed approach for intrusion detection in wireless sensor networks[END_REF] is a distributed method. The N nodes are divided into d groups of n i nodes, for i = 1, . . . , d where d i=1 n i = N . Each group i deals p i features where

where X i is a n i -by-p i submatrix of X. First, for each group, data are normalized to the range [0 , 1] and then one computes the column-centered matrix X c i . For each group i, The principal axes of X c i are extracted using a singular value decomposition:

The matrix of the principal axes is given by the matrix V i . Each group sends this matrix to its neighbor group according to a routing Finally, the principal axes of X are W = i V i . The PCADID operates through a "divide-to-conquer" scheme in order to reduce the computational complexity of the eigen-decomposition problem. However, it still requires the eigen-decomposition of several matrices, making it inappropriate for large scale networks. Moreover, our approach outperforms the PCADID technique as shown with experimental results in the following section; see for instance Fig. 2.

VII. EXPERIMENTAL RESULTS

In this section, the performance of the proposed approach is illustrated within two different contexts: dimensionality reduction of time series in wireless sensor networks and image processing. The performance is measured in terms of the angle between the principal axis w * , obtained from the centralized strategy with an eigen-decomposition of the covariance matrix, and the estimate w * ,l at node l, namely

In order to provide a fair comparative study, we use the same initial random estimate for all strategies.

A. Time series measurements in a WSN

In this section, we consider the problem of tracking a gas spread using a wireless sensor network (WSN) [START_REF] Akyildiz | A survey on sensor networks[END_REF]. The sensor k at a position denoted by z k ∈ Z measures, at time θ, a gas quantity denoted by x k,θ . The region under scrutiny Z = [-0.5, 0.5] × [-0.5, 0.5] is a two-dimensional unit-area. We aim to reduce the order of the time series of the measurements. The gas diffusion within this region is governed by the following differential equation:

where G(z, θ) is the density of gas depending on the position z and time θ, ∇ 2 z is the Laplace operator, c is the conductivity of the medium, and Q(z, θ) corresponds to the added quantity of gas. A gas source placed at the origin is activated from θ = 1s to θ = 15s. We use N = 100 sensors uniformly deployed in the region Z, each acquiring a time series of 15 measurements, between θ = 1s and θ = 15s.

We consider a predetermined range of communication in the WSN: two nodes are considered as being connected when they are less than s units of distance apart, that is

In our experiments, we set this threshold to s = 0.38. On the other hand, throughout our simulations, we have noted that the value of the second eigenvector λ 2 is generally very small (λ 2 ≈ 0). The reason behind this is that the rate of the measurement increase of each sensor, between two consecutive moments, is approximately the same for all sensors. This behavior is inherent to the nature of gas propagation in the space. Therefore, the second principal axis is difficult to be determined in this case. As for the choice of the stepsize parameters η k,t , we use a cross validation approach: we vary the stepsize according to three series of values, the first being from 0.1, 0.01, . . . , 0.1 • 10 -12 , the second from 0.25, 0.025, . . . , 0.25 • 10 -12 , and the third from 0.5, 0.05, . . . , 0.5 • 10 -12 . Then, we choose the step size which leads to the best result in the noncooperative strategy. Note that the noncooperative strategy is used at this phase because of its smaller computational cost, compared to the other strategies. We consider a stepsize η 1 = 0.0025 for the first principal axis. The same stepsize value is taken for the cooperative strategies.For the information theory, we use η 1 = 0.025, for OJAn algorithm η 1 = 0.005, for LUO algorithm η 1 = 0.0025 for both cooperative and noncooperative strategies.

Fig. 2 and Fig. 3 show the convergence of each strategy for Oja's rule, the information theory, the Rayleigh quotient, the OJAn, and the LUO algorithms respectively, in terms of the angle Θ 1 for the first axis, averaged over the N nodes. The confronted strategies are the noncooperative strategy, the combine-then-adapt (CTA) and adapt-then-combine (ATC) diffusion strategies, and the consensus implementation. For the diffusion strategies, we take the simplified case where no information exchange is performed between nodes, except for the combination step; in other words, we use c kk = 1, and c lk = 0 for each k = l, with l, k = 1, . . . , 100. Results are shown in terms of the angle averaged over all the nodes. Moreover, Fig. 2 shows that our strategies outperforms the PCADID. These learning curves show that the noncooperative strategy is outperformed by any diffusion strategy. The analysis of the diffusion strategies shows that the adapt-then-combine strategy performs slightly better that the combine-then-adapt strategy. The overall efficiency of the diffusion strategies is shown in terms of stability. On the other hand, considering the rate of convergence, the figures show that the LUO algorithm has the fastest convergence, followed consecutively by OJAn, Oja, RQ and IT algorithms. Oja and OJAn have the same behavior in terms of speed of convergence, RQ and LUO have a conflicting behavior as expected from TABLE I. Another