Paul Honeine

Mehdi Essoloh

Hichem Cédric Richard

Hichem Snoussi

Distributed regression in

come

I. INTRODUCTION

Wireless sensor networks involve large numbers of deployed tiny radio-equipped sensors, and provide an inexpensive way to monitor physical phenomena, such as a temperature field. With relatively inexpensive wireless devices, each device has a limited amount of memory, reduced processing capabilities, limited power resources, and low communication capacities. In order to carry out good coverage of the region under scrutiny, sensors must be deployed densely, resulting in highly redundant spatial information. Many researchers in the sensor network community exploit this redundancy in order to reduce the complexity of the resulting model, see for instance [START_REF] Vuran | Spatial correlation-based collaborative medium access control in wireless sensor networks[END_REF] and references therein. Guestrin et al. consider in [START_REF] Guestrin | Distributed regression: an efficient framework for modeling sensor network data[END_REF] the spatial and temporal redundancy, where the latter is handled by fitting a three-degree polynomial to the sensed data. While this is a model-based technique, model-independent approaches received considerable attention recently. Rabbat et al. applied in [START_REF] Rabbat | Distributed optimization in sensor networks[END_REF] an incremental subgradient optimization technique for parameter estimation, yielding an efficient scheme for sensor network in terms of the energy-accuracy tradeoff. As pointed out by Predd et al. in [START_REF] Predd | Distributed learning in wireless sensor networks[END_REF], such a technique is however inappropriate for functional estimation.

Recently, increasing research attention has been directed towards kernel methods for pattern recognition, in both unsupervised and supervised learning for classification and regression problems. Based on the concept of reproducing kernels initially introduced by Aronszajn in the 50's, these methods have gained popularity with the prelude of support vector machines and the statistical learning theory [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. The literature on wireless sensor networks does not escape from the proliferation of kernel machines and its attractiveness, as studied for problems such as localization [START_REF] Nguyen | A kernel-based learning approach to ad hoc sensor network localization[END_REF], detection [START_REF] Nguyen | Nonparametric decentralized detection using kernel methods[END_REF], estimation [START_REF] Snoussi | Distributed bayesian fault diagnosis in collaborative wireless sensor networks[END_REF], and regression [START_REF] Predd | Distributed kernel regression: An algorithm for training collaboratively[END_REF]. Predd et al. propose in [START_REF] Predd | Distributed kernel regression: An algorithm for training collaboratively[END_REF] to assess distributively the regression problem by solving it locally on each sensor device, which gets information from its neighbors. While this technique is computationally efficient, it has some drawbacks. On the one hand, sensor devices must be set in broadcast mode in order to communicate with their neighbors, and thus require more power for lateral communication as compared to a sensor-to-sensor communication. On the other hand, each sensor needs to solve the optimization problem, resulting in a on-sensor matrix inversion that leads to considerable computational burden. This technique suffers from these disadvantages as long as sensors are densely deployed.

Classical kernel machines are inappropriate for regression in the context of sensor network, mainly for one major drawback, illustrated by the Representer Theorem [START_REF] Kimeldorf | Some results on tchebycheffian spline functions[END_REF], [START_REF] Schölkopf | A generalized representer theorem[END_REF]: The order of the resulting model is equal to the number of available observations. Therefore, the model order scales badly with the number of deployed sensors. In order to overcome this difficulty, we consider in this paper a reduced-order approach. Controlling the complexity of the model is widely used not only in kernel machines [START_REF] Burges | Simplified support vector decision rules[END_REF], but also in sparse representation literature [START_REF] Tropp | Greed is good: algorithmic results for sparse approximation[END_REF], [START_REF] Donoho | Compressed sensing[END_REF]. It turns out that this is equivalent to selecting a small set of sensor devices that are well spread in the network. In the proposed scheme, each device updates the model from its measurement, and increments its order if necessary, then transmits the model parameters to the next device, and so on. Both this sensor-to-sensor communication scheme which reduces the overall energy consumption, and the gradient-based learning technique that we derive in this paper, allow us to overcome the disadvantages of the method proposed by Predd et al.

This paper is organized as follows. In the next section, functional learning with classical kernel machines is concisely introduced, and we illustrate via the Representer Theorem their inappropriateness for learning in wireless sensor network. A solution to this drawback is proposed by controlling the model order, as derived in Section III. Section IV is devoted to the proposed algorithm based on a gradient descent scheme, and implementation issues regarding wireless sensor networks are studied. The efficiency of the proposed approach is illustrated in Section V, with an application to estimating a spatial temperature distribution.

II. FUNCTIONAL LEARNING IN SENSOR NETWORKS

In a conventional regression problem, one seeks a function that links the best the input space to the observation domain. By designating the former by X and the latter by D, we learn this function ψ * (•) from available data, (x i , d i ) ∈ X × D for i = 1, . . . , n. The optimization problem is given by minimizing the mean-square-error between the model output ψ(x i ) and the desired output d i , with

ψ * (•) = arg min ψ 1 n n i=1 |d i -ψ(x i )| 2 .
It is well-known that one needs to constrain this optimization problem into a functional space of smooth functions. This can be done by combining the reproducing kernel Hilbert space (RKHS) with Tikhonov regularization, yielding the optimization problem

ψ * (•) = arg min ψ∈H 1 n n i=1 |d i -ψ(x i )| 2 + η ψ 2 H , (1)
where the parameter η controls the tradeoff between smoothness and fitting the data. In this expression, H is the RKHS of a given reproducing kernel κ(•, •). This means that each function of H is evaluated at any x ∈ X with ψ(x) = ψ(•), κ(•, x) H , where •, • H is the inner product in H. Initially introduced in [START_REF] Kimeldorf | Some results on tchebycheffian spline functions[END_REF] and generalized more recently in [START_REF] Schölkopf | A generalized representer theorem[END_REF] for the wide class of kernel machines, the Representer Theorem states that the solution of the regularized optimization problem (1) is of the form

ψ * (•) = n k=1 α k κ(x k , •). (2)
In order to determine the model, completely identified by its coefficients α k , we inject this expression back into (1), yielding the so-called dual optimization problem

α * = arg min α d -Kα 2 + η α ⊤ Kα, (3)
where d and α are n-by-1 column vectors whose i-th entries are d i and α i , respectively, and K is the n-by-n matrix whose (i, j)-th entry is κ(x i , x j ). This is a linear optimization problem, and its solution is given by

α * = K ⊤ K + ηK -1 K ⊤ d. (4)
Consider the regression problem in sensor networks where, for instance, we seek to estimate a temperature field. The input space corresponds to the location in the region under scrutiny, X ⊂ IR 2 in a plane, and the output is the temperature measure D ⊂ IR. Sensors measure temperature at their respective locations, and therefore construct the training set. In a distributive scheme, each sensor updates the coefficient vector with available information and transmits it to the next sensor.

However, solving the optimization problem described above can become cumbersome as illustrated by the matrix inversion in (4), whose computational complexity is O(n 3 ), and therefore scales badly with the number of sensors in the network. While we derive this result for the particular case of the mean-square-error as a cost functional, such drawback is common to all classical kernel machines as illustrated by the Representer Theorem, with a model completely determined by the n coefficients, α k , as well as the n coordinates x k representing sensor locations. In order to propose a distributive technique to solve this problem, Predd et al. propose in [START_REF] Predd | Distributed kernel regression: An algorithm for training collaboratively[END_REF] to solve it locally by every sensor, with location and measurement of each neighboring sensor taking part in the calculations. This requires to convey this information and to invert a matrix similar to the one in (4) whose dimensions correspond to the number of available neighbors. Such technique is however unadapted in practice as sensors are densely deployed, with heavier communication burden and higher neighborhood concentration.

III. A REDUCED-ORDER MODEL FOR SENSOR NETWORKS

In order to overcome the difficulties invoked in the previous section, we propose a reduced-order model constructed from the optimal model (2) by considering only a small number of kernel functions in the expansion. Let m be that number, the reduced-order model is defined by

ψ(•) = m k=1 α k κ(x ω k , •), (5)
where the m coordinates x ω k are selected from the locations of all the sensors, and thus ω k ∈ {1, . . . , n}. Without going into details, we note that by injecting (5) back into the cost functional (1) of the regularized optimization problem, we obtain a solution similar to [START_REF] Predd | Distributed learning in wireless sensor networks[END_REF], where one needs to inverse an m-by-m matrix with entries of the form κ(x ωi , x ωj ). In Section IV, we derive a receive-update-transmit scheme for each sensor in order to learn the model by minimizing the error on its sensed information. But before, we study the problem of selecting these kernel functions, or equivalently selecting the m sensors with locations given by x ω1 , . . . , x ωm . The problem of sparse representations spans many scientific and research area, and more specifically in signal processing [START_REF] Tropp | Greed is good: algorithmic results for sparse approximation[END_REF], [START_REF] Donoho | Compressed sensing[END_REF] and pattern recognition with kernel machines [START_REF] Hoegaerts | Eigenspace methods and subset selection in kernel based learning[END_REF], [START_REF] Csató | Sparse representation for gaussian process models[END_REF]. Sparse techniques with ℓ 1 -based penalization are too computational expensive for wireless sensor networks. A more appropriate approach would be the one proposed in [START_REF] Baudat | Kernel-based methods and function approximation[END_REF], [START_REF] Csató | Sparse representation for gaussian process models[END_REF], which translates in sensor networks as follows: the model is constructed in a simple walk through the network, where each sensor discards its kernel function from the model, and thus leaves its order unchanged, if the kernel function can be well approximated by the model; otherwise the order is incremented by adding the kernel function to the model. In other words, at sensor i, if the quantity min β1,...,βm

κ(x i , •) - m k=1 β k κ(x ω k , •) 2 H
is smaller than a given threshold, then the kernel function κ(x i , •) is not included into the model; otherwise, it is added to the model. This criterion however requires the inversion of an m-by-m matrix, and therefore demands high precision and computational resources from the on-sensor microprocessor for each sensor.

We propose in this article to reduce further the computational burden of evaluating the criterion for determining the relevance of a kernel function with respect to the model. For this, we take advantage of recent developments investigated by two of the authors, on the coherence criterion defined as follows: The kernel function κ(x i , •) is included to the morder model if

max k=1,...,m |κ(x i , x ω k )| ≤ ν, (6)
where ν is a given threshold. This means that the kernel functions of the resulting model have a bounded cross-correlation, since from ( 6) we have

max k,l=1,...,m | κ(x ω k , •), κ(x ω l , •) | ≤ ν.
Without going into details, properties of this criterion and connections to other sparsification criteria are studied in our recent papers [START_REF] Honeine | On-line nonlinear sparse approximation of functions[END_REF], [START_REF] Richard | Online prediction of time series data with kernels[END_REF] for prediction in time series data. From this criterion, we derive a natural one for learning in wireless sensor networks.

In the particular case of the widely used radial kernels which can be expressed in terms of κ( x ix j ), the criterion (6) may be simplified further. This is because κ(x i , x ω k ) can be written in terms of x ix ω k . To illustrate this property without loss of generality, we consider the Gaussian kernel defined by κ(x i , x j ) = e -xi-xj 2 /2β 2 0 where β 0 is a tunable parameter. By substituting this definition into the criterion expression above, we obtain the expression min k=1,...,m

x i -x ω k 2 > 2β 2 0 ln(1/ν).
In other words, sensor i increments the model order by including its kernel function to the model, if it does not belong to the neighborhood of sensors previously in the model. The notion of neighborhood here is defined by the right-hand-side of the above expression. By setting ν 2 0 = 2β 2 0 ln(1/ν), we get the so-called neighborhood criterion min k=1,...,m

x i -x ω k > ν 0 . (7)
Thus, it turns out that the criterion (6) based on functional approximation in the RKHS is natural for the topology of sensor networks, as it is equivalent to the neighborhood criterion defined in [START_REF] Nguyen | Nonparametric decentralized detection using kernel methods[END_REF]. Next we derive an algorithm to learn such reduced-order model, by considering the criterion [START_REF] Nguyen | Nonparametric decentralized detection using kernel methods[END_REF].

IV. DISTRIBUTED REDUCED-ORDER MODEL ALGORITHM

A large class of algorithms can be proposed for solving the optimization problem, mainly with techniques studied in [START_REF] Honeine | On-line nonlinear sparse approximation of functions[END_REF], [START_REF] Richard | Online prediction of time series data with kernels[END_REF] and adapted here to the neighborhood criterion. In what follows, we derive a low-computational demanding algorithm, based on a simple gradient descent approach, and study implementation issues in wireless sensor networks.

A. The algorithm

For this purpose, each sensor i updates the coefficient vector α i to be as close as possible to the previous one α i-1 by imposing a null error on approximating the measurement. This is equivalent to solving the following constraint optimization problem :

min αi α i-1 -α i 2 (8) subject to κ ⊤ i α i = d i , (9)
with κ i an m-by-1 column vector whose k-th entry is κ(x ω k , x i ), where the constraint is obtained by applying (5) to information relative to sensor i. Solving this problem at each sensor depends on the neighborhood criterion, and thus two cases are possible, depending if ( 7) is satisfied or not.

Case 1. min k=1,...,m x ix ωk < ν 0 :

This is the case when sensor i is in the neighborhood of any of the m previously selected sensors ω 1 , . . . , ω m . Therefore its kernel function may be well approximated by a linear combination of the model kernel functions. To solve the constraint optimization problem ( 8)-( 9), we consider minimizing the corresponding Lagrangian given by

α i-1 -α 2 + λ(d i -κ ⊤ i α),
where λ is the Lagrangian multiplier. By setting to zero the derivatives of the above cost function with respect to α and λ, we get the following conditions on α i to verify:

2 (α i -α i-1 ) ⊤ = λ κ ⊤ i and κ ⊤ i α i = d i . Since κ ⊤ i κ i is nonzero, these equations lead to λ = 2 (κ ⊤ i κ i ) -1 (d i -κ ⊤ i α i-1 )
, where d i -κ ⊤ i α i-1 is the a priori error at sensor i. By substituting into the above expression, we obtain the update equation

α i = α i-1 + ρ κ i 2 κ i d i -κ ⊤ i α i-1 . ( 10
)
where we have introduced the step-size control parameter ρ, as preconised in conventional adaptive filtering techniques.

Case 2. min k=1,...,m x ix ωk > ν 0 :

In this case, the sensor i is not covered by the m sensors defining the kernel functions of the model. Therefore, we increment the model order by including the corresponding kernel function to it. To accommodate the new element in m + 1 vector α i , we modify problem ( 8)-( 9) as

min α α [1,...,m] -α i-1 2 + α 2 m+1 subject to d i = κ ⊤ i α,
where α [1,...,m] denotes the first m elements of the vector, and κ i has been increased by one entry. Considerations similar to those made in Case 1 lead to the updating rule

α i = α i-1 0 + ρ κ i 2 κ i d i -κ ⊤ i α i-1 0 , (11)
with κ i the (m + 1)-by-1 column vector whose k-th entry is κ(x ω k , x i ), and ω m+1 = i.

[x 1 ] [x 1 ] [x 1 ] [α 1,1 ] [α 1,2 ] [α 1,3 ] [x 1 x 4 ] [x 1 x 4 ] [α 1,4 α 2,4 ] [α 1,5 α 2,5 ]
Fig. 1. A schematic representation illustrating the proposed approach.

B. Analysis of the algorithm

We illustrate the proposed approach in Fig. 1, where we present five sensors with links representing the walk of the information throughout the network. Next we denote by α k,i the k-th entry of the coefficient vector α i updated by sensor i. Initializing at Sensor 1 yields ω 1 = 1 and m = 1. Its location x 1 and the estimated coefficient α 1,1 are transmitted to the next sensor in the network. Sensor 2 belongs to the neighborhood of Sensor 1, since x 2x ω1 is less than the given neighborhood threshold. Thus it updates the coefficient into α 1,2 with the rule (10) by leaving the model order unchanged. It transmits α 1,2 and x 1 to the next sensor, Sensor 3, with updating similar to Sensor 2 since it also belongs to the neighborhood of Sensor 1. Not in that neighborhood, Sensor 4 increments the model order by adding its kernel function to the model and a new weighting coefficient. Updating with [START_REF] Schölkopf | A generalized representer theorem[END_REF], the resulting coefficients as well as the locations [x 1 x 4 ] are transmitted to the next sensor, and so on.

In designing algorithms for wireless sensor network applications, several implementation issues must be taken into consideration, mainly algorithmic complexity and energy consumption, the latter being proportional to the size of transmitted messages. The proposed approach answers this question, on both aspects. On the one hand, the gradient technique has low computational complexity and minor memory requirements, with no matrix inversion as opposed to Predd's technique [START_REF] Predd | Distributed learning in wireless sensor networks[END_REF], and thus can be efficiently implemented on tiny low-cost processors. On the other hand, the messages communicated between sensors are constituted only of the m coordinates and m coefficients taking part in the reduced-order model, as opposed to the full-order model where all the coordinates and coefficient values must be conveyed between sensors.

It is worth noting that for small values of the neighborhood threshold, the influence region of each sensor is reduced, resulting into a model with a larger order, therefore closer to the optimal full-order model given in [START_REF] Guestrin | Distributed regression: an efficient framework for modeling sensor network data[END_REF]. However, the price to pay is an increasing energy consumption, mainly due to the communication for conveying the parameters of the model. This is the classical accuracy-energy dilemma in wireless sensor networks, with the neighborhood threshold determines the tradeoff between the optimal approximation and the reduction in the model order.

V. SIMULATIONS

In order to illustrate the proposed approach, we consider estimating the spatial temperature distribution, with sensors densely and randomly spread over the region under scrutiny. Due to lack of experimental data from densely deployed sensor networks, we consider simulation results 1 obtained from classical partial differential equations. We consider three heat sources randomly set inside a square region, and estimate the temperature distribution from a single measurement available on each of the n = 100 sensors randomly spread.

As the configuration settings determines the neighborhood threshold to be used, we set in what follows ν 0 = 0.3, which yields a 15-order model, corresponding to a 85% of reduction in the model order compared to the full-order one. In what follows, the Gaussian kernel is used. To determine the influence of both its bandwidth and the step-size parameter, 1 A detailed description on the experimental configuration and the dataset is available in http://www.ulb.ac.be/di/mlg/sensorNet/modeling.html. we define the root-mean-square (rms) error by the square root of

1 n n i=1 (d i -κ ⊤ i α i-1 ) 2
, where d i -κ ⊤ i α i-1 = ǫ i is the a priori error of the sensor i. These results are sketched in Fig. 3 for different values of the step-size parameter. As shown, the optimal bandwidth is given by β 0 = 0.2, independently of the chosen step-size parameter and with almost the same rms error, except for ρ = 0.1. Such exception is mainly due to low convergence resulting from small values of ρ. In order to circumvent such problem, one may increase the number of sensors in the network, or more likely consider several cycles through the network, as proposed in [START_REF] Rabbat | Distributed optimization in sensor networks[END_REF]. However, this case is beyond the scope of this work, as we stick to the convergence case given by any of the other studied values of ρ.

To illustrate this convergence behavior, we plot in Fig. 4 the evolution of the (absolute) values of the a priori error as we visit each of the sensors, for cases ρ = 0.1 and ρ = 0.7. For the latter case, the convergence occurs within the first couple of sensors, as shown in the figure. The resulting temperature distribution is illustrated in Fig. 5, where we distinguish sensors contributing to the incrementation of model order. As expected, these sensors span the region under scrutiny.

VI. CONCLUSION

In this work, we showed that regression in wireless sensor networks can take advantage of recent developments in sparse representations with kernel machines. In order to derive a reduced-order model, we developed a criterion natural to the topology of the deployed sensors. A gradient-based algorithm is proposed for this purpose, and implementation issues are discussed. Simulation results show good convergence of the proposed approach.

Current and future works concentrate on expanding this approach for monitoring the evolution of the temperature with time. One way to achieve this is to use a reproducing kernel that takes into account both the spatial and the temporal information, for instance by combining kernels defined on each of those domains.

Fig. 2 .

 2 Fig. 2. Map of the sensors, randomly (and uniformly) distributed over the region under scrutiny.

Fig. 3 .

 3 Fig. 3. Evolution of the error as a function of the bandwidth of the Gaussian kernel, for different values of the step-size parameter.

Fig. 4 .

 4 Fig. 4. Convergence of the algorithm as we walk through the 100 sensors of the network.

Fig. 5 .

 5 Fig. 5. Resulting temperature map obtained from the Gaussian kernel with (ρ = 0.7, β 0 = 0.2). The three heat sources are represented by big yellowdiscs •, sensors by blue-dots •, and incrementing-order ones by red-discs • •.


Paul Honeine
email: paul.honeine@utt.fr

Cédric Richard
email: cedric.richard@utt.fr

José Carlos

M Bermudez
email: j.bermudez@ieee.org

Hichem Snoussi
email: hichem.snoussi@utt.fr

Distributed prediction of time series data with kernels and adaptive filtering techniques in sensor networks

Wireless sensor networks are becoming versatile tools for learning a physical phenomenon, monitoring its variations and predicting its evolution. They rely on low-cost tiny devices which are deployed in the region under scrutiny and collaborate with each other. Limited computation and communication resources require special care in designing distributed prediction algorithms for sensor networks. In this communication, we propose a nonlinear prediction technique that takes advantage of recent developments in kernel machines and adaptive filtering for online nonlinear functional learning. Conventional methods, however, are inappropriate for large-scale sensor networks, as the resulting model corresponds to the number of deployed sensors. To circumvent these drawbacks, we consider a distributed control of the model order. The model parameters are transmitted from sensor to sensor and updated by each sensor based the measurement information. The model order is incremented whenever this increment is relevant compared to a reduced-order model. The proposed approach is naturally adapted for predicting a varying phenomenon, as model order increases are governed by the novelty of the new observation at each sensor node. We illustrate the applicability of the proposed technique by some simulations on establishing the temperature map in an region heated by sources.

I. INTRODUCTION

Wireless sensor networks (WSN) are becoming a multidisciplinary research area attracting researchers from microelectronics, communication and signal processing communities, only to name a few. Mainly motivated by military applications, they are appropriate for a broad class of applications, including monitoring physical or environmental conditions, such as temperature, humidity or pollutants for supervising an agricultural field at different locations. They rely on tiny devices with sensing capabilities, communicating wirelessly with each other. In order to cover a large spatial region of interest, a big number of these devices (simply denoted as sensors hereafter) are deployed randomly, restricting the manufacturing cost of each one. This leads to resource-constraint sensors, limiting the complexity of the algorithm run on each one.

Sensors are mainly deployed to get some measurements at known (random) locations 1 in the region under scrutiny. Since communication between sensors is possible, at least between neighbors, they may collaborate with each other to learn the overall physical phenomenon in each location in the region, from the measurements as well as the locations of the sensors. As opposed to a fusion center approach where sensors only sense and convey information, here we have smart sensors. A common learning scheme consists of updating the model, at each sensor in turn, with the newly acquired data, and forwards information to another nearby sensor.

One can take advantage of adaptive learning, as preconized by the vast literature on adaptive filtering techniques [START_REF] Haykin | Adaptive filtering theory[END_REF], [START_REF] Sayed | Fundamentals of adaptive filtering[END_REF]. However, most of these are linear signal processing techniques, inappropriate for the inherent nonlinear nature of the wide class of systems under investigation. Nonlinear techniques are becoming more widely investigated thanks to recent developments in machine learning, and more specifically with kernel machines; See for instance [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF] for a review of kernel methods. Most of these techniques are nonlinear analysis and processing techniques obtained by revisiting linear ones in the light of reproducing kernel Hilbert spaces. There is a growing interest around this concept for data mining and machine intelligence, since recent burst of Vapnik's support vector machines and statistical learning theory [START_REF] Vapnik | Statistical Learning Theory[END_REF]. While most kernel machines yield optimal solutions, they are not established in the same line as online learning, as the order of the resulting model grows along with the number of available observations.

In this paper, we propose a nonlinear kernel-based adaptive method, applied here for WSN, by bridging the gap between the nonlinear kernel machines and adaptive filtering techniques. In [START_REF] Honeine | On-line nonlinear sparse approximation of functions[END_REF], [START_REF] Richard | Online prediction of time series data with kernels[END_REF], we have proposed to control the order of the resulting model, with the pre-processing kernelcoherence criterion, and revisited in [START_REF] Honeine | Distributed regression in sensor networks with a reduced-order kernel model[END_REF] for learning with WSN. Nevertheless, such criterion is used prior to any learning scheme, and is independent of the prediction error. Here, we revisit this criterion to propose a post-processing criterion, the functional-coherence criterion, which is applied to the learnt function, thus taking into account the prediction error. The rest of this paper is organized as follows. In section II, learning in a RKHS is briefly introduced, with focus on functional leaning. In section III, we derive techniques to control the order. Then the proposed method is studied in section IV. Simulation results are reported in section V, and finally, conclusion remarks are summarized in section VI.

II. FUNCTIONAL LEARNING WITH KERNEL MACHINES

In order to understand the nature of functional learning with kernel machines, it will be useful first to take a look at a conventional least-squares problem for fitting data. Let (x 1 , d 1 ), (x 2 , d 2 ), . . . , (x n , d n ) be a set of n training pairs, where x i is an input example drawn from a vector-space X , and d i ∈ D is the desired observation or label associated to it. The goal is to learn a map from X to D given the training set, and therefore getting the appropriate d for any x. Mathematically, the mapping or system connecting the two spaces is modeled by a function, say ψ. When the labels take values in a finite set, the task is called classification or discrimination. In a more general setting 2 , when D = IR we have a regression problem where d i corresponds to the desired output of the system ψ for the input x i . The optimal function ψ * is obtained by minimizing some cost functional on the available training set, such as the mean square error, by solving the optimization problem

ψ * = arg min ψ n i=1 |ψ(x i ) -d i | 2 .
(1)

However, this is an ill-posed problem, as an infinite number of functions sets to zero this cost functional. In order to obtain a well-posed problem, one must restrict the space of candidate functions into for instance linear, quadratic, or polynomial functions or even splines. Linear functions are conceptually easy and simple to implement, as preconized by the large literature on adaptive filtering techniques [START_REF] Haykin | Adaptive filtering theory[END_REF], [START_REF] Sayed | Fundamentals of adaptive filtering[END_REF]. For this purpose, we restrict ourselves to functions ψ defined by a vector a ∈ X such that ψ(x) = a, x for any x ∈ X , where •, • is the usual inner product in a vector space. Thus, the optimization problem is reduced to determine the vector a * by solving

a * = arg min a n i=1 a, x i -d i 2 , (2)
This optimal vector can be learnt recursively in an online learning scheme as suggested with adaptive filtering techniques, where at each time-instance n we have a new observation pair (x n , d n ). Nevertheless, this optimization problem may still be ill-posed, and some structure of the candidates vectors a must be required, taking into account for instance the regularity of the system as well as the noisy measurements. Moreover, it is worth noting that, by their nature, such linear functions are not suitable to nonlinear analysis and prediction.

In this paper, we take advantage of recent developments in the area of nonlinear functional learning, with kernel machines. The space of candidate functions is constructed on the underlying structure of reproducing kernel functions, as defined next. Let (H, • , • H ) denotes a Hilbert space of realvalued functions defined on X . A reproducing kernel of

H Polynomial κ(x i , x j ) = ( x i , x j + β 2 ) p Laplace κ(x i , x j ) = exp `- x i -x j /σ 2 0 Ǵaussian κ(x i , x j ) = exp `- x i -x j 2 /2σ 2 0 TABLE I SOME REPRODUCING KERNELS, WITH TUNABLE PARAMETERS β, p, σ 0 .
is a function κ(x i , x j ) from X × X to IR that verifies both properties:

• the function κ(x i , •), belongs to H for all x i ∈ X , • ψ(x i ) = κ(x i , •), ψ H for all x i ∈ X and ψ ∈ H.
The Hilbert space H associated with this kernel is called a reproducing kernel Hilbert space (RKHS). The last property is a fundamental property of RKHS, which states that the evaluation of a function ψ ∈ H in that space at any x i ∈ X is given by its inner product with the function κ(x i , •). Moreover, one can write for the particular case ψ = κ(x j , •) for any

x j ∈ X the property κ(x i , •), κ(x j , •) H = κ(x i , x j ).
Therefore, κ(x i , x j ) corresponds to an inner product between data in a modified or transformed space. Back to the optimization problem (1) and its linear counterpart (2) applied to the RKHS, we get the following optimization problem

ψ * = arg min ψ∈H n i=1 |ψ(x i ) -d i | 2 + η ψ 2 H . (3)
In this expression, we use a Tikhonov regularization with η ∈ [0, 1], allowing a tradeoff between the well-fitness to the data (the first term in the above cost function) and the regularity of the resulting function (the second term). It is worth noting that we get a linear-function problem (2) by considering the linear reproducing kernel κ(x i , x j ) = x i , x j . Nonlinearity is achieved by substituting the inner product by a nonlinear reproducing kernel as defined in Table I, leaving the algorithm unchanged and incurring essentially the same computational cost.

III. MODEL ORDER CONTROL FOR ONLINE LEARNING

By virtue of the well known Representer Theorem [START_REF] Kimeldorf | Some results on tchebycheffian spline functions[END_REF], [START_REF] Schölkopf | A generalized representer theorem[END_REF], solutions to the regularized optimization problem (3) have the form

ψ * = n j=1 α * j κ(x j , •), (4)
where the order of the model corresponds to the number of available data. By analogy to the linear case defined in (2), x i are substituted here by their corresponding kernel functions κ(x i , •), and a * by the function ψ * defined in [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF]. Inserting this model into (3) leads to the following optimization problem:

α * = arg min α d -Kα 2 + η α ⊤ Kα,
where α and d are n-dimensional column vectors of entries α k and d k , respectively, and K is the n-by-n (Gram) matrix of entries κ(x i , x j ), for i, j = 1, . . . , n. From classical matrix algebra, we get the solution to this well known problem with

(K + ηI) α * = d,
with I the identity matrix of appropriate size. Therefore, this requires the inversion of the matrix K + ηI, which is warrantable thanks to a non-zero regularization constant η.

In an online setting, we have a new data pair (x n , d n ) at each instant n, and the model is updated recursively based on this new information. However, considering a model of the form (4), its order will steadily increasing. In order to circumvent this inconvenience, we propose to control its order. This is done by retaining only a small number of kernel functions in the expansion (4), say m, and write

ψ n = m k=1 α n,k κ(x ω k , •), (5)
where ω 1 , . . . , ω m are the indexes of the retained kernel functions, from the n available ones. The selection criterion determines if, at instant n, the kernel function κ(x n , •) should be included to the expansion, or discarded, or even removed if it already belongs to it. Therefore, the model order may increment, stay unchanged, or decrement. In [START_REF] Honeine | On-line nonlinear sparse approximation of functions[END_REF], [START_REF] Honeine | Distributed regression in sensor networks with a reduced-order kernel model[END_REF], [START_REF] Richard | Online prediction of time series data with kernels[END_REF],

we derived an easy to compute criterion applied in a preprocessing scheme, independently of the resulting optimal function as well as the prediction error. In order to circumvent these drawbacks, we investigate in this paper a new criterion applied in a post-processing scheme, after getting the optimal function and therefore is observation-aware as opposed to the former. Since this criterion follows the same spirit of the former, we begin by a brief review of the former criterion.

A. The kernel-coherence criterion, a pre-processing criterion

In this first approach, we seek a compact representation, where the current kernel function is discarded from the model expansion when it can be well approximated by the already available ones. This can be done by comparing κ(x n , •) to its projection onto the space spanned by the other m kernel functions, yielding the following optimization problem

min γ κ(x n , •) - m j=1 γ j κ(x ωj , •) 2 H < ν 2 0 , (6)
for a given threshold ν 0 . From matrix algebra, this is solved by an inversion of the m-by-m Gram matrix of entries κ(x ωi , x ωj ) for i, j = 1, . . . , m, as given for instance in [START_REF] Engel | The kernel recursive least squares algorithm[END_REF], and thus has a computational complexity of O(m 3 ). In [START_REF] Honeine | On-line nonlinear sparse approximation of functions[END_REF], [START_REF] Richard | Online prediction of time series data with kernels[END_REF], we consider the following more efficient criterion. We discard κ(x n , •) from the expansion if, for a given threshold ν 1 , we have max

j=1,...,m | κ(x n , •), κ(x ωj , •) H | κ(x n , •) H κ(x ωj , •) H > ν 1 . (7)
In other terms, this corresponds to upper-bounding the cosine of the angles between the kernel functions in the resulting expansion. For a small threshold value, we get a set of incoherent kernel functions, and orthogonality when ν 1 = 0; we call this criterion the kernel-coherence criterion. This criterion is very efficient from a computational point of view, since the numerator of ( 7) is given by |κ(x n , x ωj )| and its denominator is κ(x n , x n )κ(x ωj , x ωj ).

In [START_REF] Richard | Online prediction of time series data with kernels[END_REF], we derive the connection between both criteria (6) and [START_REF] Richard | Online prediction of time series data with kernels[END_REF]. These criteria consider the most compact representation (9) that approximates well the optimal model (4). However, they do not take into account the observation d n associated to κ(x n , •), neither the prediction error. This is mainly due to the pre-processing nature of the approach, as selection occurs before any computation of the resulting model, with and without the kernel function κ(x n , •).

B. The functional-coherence criterion, a post-processing criterion

The kernel-coherence criterion mentioned above, as defined in [START_REF] Richard | Online prediction of time series data with kernels[END_REF], determines the relevance of the information in κ(x n , •), with respect to all previous kernel functions in the expansion. We propose to revisit this principle, and apply it to the learnt function. For now, consider that the function ψ n is computed with an expansion including κ(x n , •). We propose to reduce the model order by discarding this kernel function if ψ n is not very relevant, compared to previous computed functions, with max k=1,...,n-1

ψ n , ψ k H ψ n H ψ k H > ν ′ 0 ,
where ψ k is the optimal function obtained at instant k. This criterion may be considered as an extension of [START_REF] Richard | Online prediction of time series data with kernels[END_REF], and in the same line of upper-bounding the cosine of the angles between functions in order to span the largest subspace in the RKHS.

As opposed to [START_REF] Richard | Online prediction of time series data with kernels[END_REF], this criterion considers the whole learnt function, including its prediction error, and not only its kernel expansion. While this requires keeping all previously learnt functions, we derive a sufficient condition by considering the criterion

ψ n , ψ ⊥ n H ψ n H ψ ⊥ n H > ν ′ 1 ,
for a given threshold ν ′ 1 ∈ [0, 1], and where ψ ⊥ n is the orthogonal projection of ψ n onto the space spanned by the n -1 functions. Moreover, we don't need to keep all these functions in memory, since this space is also spanned by the retained kernel functions. In the next section, we propose a learning algorithm, where this criterion is applied at the end in order to determine which of ψ n or ψ ⊥ n must be considered as the resulting optimal function, based on this measure of coherence between them.

IV. DISTRIBUTED FUNCTIONAL LEARNING IN WSN

Consider the problem of functional learning in wireless sensor networks (WSN), such as estimating a temperature field for instance, where the input of the function is the coordinates in the region under scrutiny, and its output is estimated temperature. We learn the optimal function from sensors deployed in this region, giving us measurements at known locations. Since most considered techniques are iterative, this paper is not an exception, more than one pass through the network may be required in order to converge to the optimal solution, as well as in case of tracking the evolution of the system.

In [START_REF] Honeine | Distributed regression in sensor networks with a reduced-order kernel model[END_REF], we derived a distributed regression technique based on the kernel-coherence criterion. However, this criterion depends only on the location of the sensors, x n , thus the model is fixed for a given deployment for a the large class of motionless-sensor problem (at least after the first pass through the network), and only the weighting coefficients are updated. The functional-coherence criterion allows us to overcome this drawback, since it looks at the functions where the prediction error is considered in the iterative updating scheme. This approach is illustrated in this paper, with a kernel-based normalized least-mean-squares (NLMS) technique derived next.

Consider the optimization problem (3), where ψ is substituted by the reduced-order model [START_REF] Kimeldorf | Some results on tchebycheffian spline functions[END_REF]. Each sensor, in turn, updates the model ψ n with the newly acquired data d n and forward information to another nearby sensor, and so on. The following actions are performed:

A. An updating parameter stage, with or without changing the model order depending on whether κ(x n , •) does not belong to the the model expansion, or it already belongs. The model resulting from this stage, ψ n incorporates this kernel function. Let m + 1 be its order. B. A projection stage, where we compute ψ ⊤ n , the projection of ψ n onto the space with one removed kernel function, κ(x n , •). The resulting model, ψ ⊥ n does not include the latter. Its order is m. C. A model order reduction stage, by considering as a final model ψ ⊤ n rather than ψ n , if

ψ n , ψ ⊥ n H ψ n H ψ ⊥ n H > ν ′ 1 .
Next, we develop each of these stages, while the first stage is decomposed into two cases, depending if the current kernel function already belongs to the model, or not. For the sake of simplicity, we begin by developing the algorithm for η = 0 in both cases, and introduce the regularization parameter at the end of each.

A1. if κ(x n , •) already belongs to the model

Without changing the model order, defined by ψ n = m+1 k=1 α n,k κ(x ω k , •), we consider updating the coefficients, α n,k , by solving the following minimum-deviation problem

min α α -α n-1 2 subject to α ⊤ κ n = d n , where κ n = κ(x ω1 , x n ) • • • κ(x ωm+1 , x n )
⊤ , and the resulting α n is a column vector of entries α n,k for k = 1, . . . , m+1. This constraint is equivalent to

ψ n (x n ) = d n .
To solve the constraint optimization problem, we write the Lagrangian

α -α n-1 2 + λ(d n -α ⊤ κ n )
where λ denotes the Lagrangian multiplier. By taking to zero its derivatives with respect to α and λ, we get the two conditions on α n :

α n -α n-1 = 1 2 λ κ ⊤ n and α ⊤ n κ n = d n From these conditions on α n-1 , we get λ = 2(κ ⊤ n κ n ) -1 (d n - α ⊤ n-1 κ n )
, where κ ⊤ n κ n is assumed non-singular. By injecting this expression in the first condition, we get the recursive update

α n = α n-1 + ρ η + κ n 2 κ n (d n -α ⊤ n-1 κ n ), (8)
where ρ is step-size control parameter, preconised in the adaptive filtering literature, and where we reintroduce the regularization parameter η. In this expression, d nα ⊤ n-1 κ n corresponds to the prediction error resulting from the previous model.

A2. if κ(x n , •) does not belong to the model

In this case, we increment the model order by including this kernel function to the expansion, yielding the model

ψ n = m k=1 α n,k κ(x ω k , •) + α n,m+1 κ(x n , •). (9)
With one more coefficient and kernel function, the minimumdeviation problem becomes

min α α [1:m] -α n-1 2 + α 2 m+1 subject to α ⊤ κ n = d n ,
where α [1:m] denotes the first m entries of vector α

= [α 1 • • • α m α m+1 ] ⊤ , and κ n = [κ(x ω1 , x n ) • • • κ(x ωm , x n ) κ(x n , x n )] ⊤ .
Similarly to A1, we get the recursive updating

α n = α n-1 0 + ρ η + κ n 2 κ n d n -κ ⊤ n α n-1 0 . (10)

B. Projection stage

In order to simplify the notations, let κ(x ωm+1 , •) be the current kernel function, i.e. n = ω m+1 Let ψ ⊥ n be the orthogonal projection of ψ n onto the space spanned by the other m kernel functions, obtained by solving κ(x ωj , •), ψ n -ψ ⊥ n H = 0 for all j = 1, . . . , m.

Since we have ψ ⊥ n = m j=1 β n,j κ(x ωj , •), this is done by solving the m-by-m linear system

K ω β n = ψ ω , (11)
where β n and ψ ω are m-element column vectors of entries β n,j and ψ n (x ωj ) for j = 1, . . . , m, respectively, and K ω is the m-by-m matrix of entries κ(x ωi , x ωj ) for i, j = 1, . . . , m.

It is worth noting that we can make use of a rank one update for the matrix inversion process.

C. Model order reduction stage

In both cases, A1. and A2., the resulting model includes the kernel function κ(x n , •). Now we are in a position to reduce the model order, by using ψ otherwise, we keep ψ n as the final model. Let α n be the column vector of the weighting coefficients, obtained from ( 8) or [START_REF] Schölkopf | A generalized representer theorem[END_REF]. Then, this criterion is equivalent to

⊥ n instead of ψ n if ψ n , ψ ⊥ n H ψ n H ψ ⊥ n H > ν ′ 1 ;
(α ⊤ n ψ ω ) 2 α ⊤ n K ω α n ψ ⊤ ω K -1 ω ψ ω > ν ′2 1 ,
where relation ( 11) is used.

V. SIMULATIONS

In order to illustrate our approach, we consider the problem of monitoring diffusion phenomena governed by the generic partial differential equation

∂C(x, t) ∂t -∇(D ∇C(x, t)) = Q(x, t)
with C the concentration as a function of location and time, Q a volume source, which may depend on both location and time, and the diffusion parameter D (scalar or matrix), depending on the medium. When D is a constant, the diffusion is isotropic, and when D varies with direction, it is called anisotropic. In this paper, we consider an anisotropic medium 3 , with several concrete examples such as studying the behavior of a pollutant concentration in the air in a windy environment, or monitoring a contaminating agent in water resources where the water flow takes part in its spread, in a given direction.

For this purpose, we consider 100 sensors in a fixed 2D location, deployed randomly in the 1.5-by-1.5 square region under scrutiny. We study the evolution of the system from t = 1 to t = 20, where at each instance we have one new measurement per sensor. From results on preliminary experiments done on instant t = 10, with 10 passes, we set η = 0, ρ = 0.99 and ν ′ 1 = 0.99. The commonly investigated Gaussian kernel is used, κ(x i , x j ) = exp(x ix j 2 /2σ 2 ), with its bandwidth set to σ = 0.2.

Figure 1 illustrates the evolution of the resulting model, at different time instances. We see that at time instant t = 1, we have a concentration within a given region. While the 3 For simulations, we use the PDE toolbox of Matlab, with D a function increasing with the x-coordinate, and Q a time-increasing source. diffusion operates mainly in the x-direction, our model follows this anisotropic diffusion. This is illustrated by the contour lines in Figure 1, as well as with the (big-red) selected sensors defining the model, spanning the region under evolution VI. CONCLUSION This paper introduces a new criterion to control the model order of kernel machines. Compared to previously proposed criterion which are applied prior to any learning task and independent of the observations, this functional-coherence criterion operates on the learnt function. It can be regarded as a measure of coherence between the learnt function and it projection onto a smaller-dimension space, which determines if we keep the function or its projection as the final result.

Fig. 1 .

 1 Fig. 1. The evolution of the model at t = 1, t = 2, and t = 20, with contour lines. The blue dots and the big-red dots represent the sensors and those contributing to the model, respectively.

Even when deployed randomly, we assume that each sensor is locationaware, with coordinates obtained by a self-localization pre-processing technique, as illustrated in[START_REF] Honeine | Localization in sensor networks -a matrix regression approach[END_REF] and references therein. Moreover, most applications involve motionless devices.

An even more general setting consists of D = IR d or even D = l C d , however such straightforward generalization is beyond the scope of this communication.


Mehdi Essoloh
email: mehdi.essoloh@utt.fr

Cédric Richard

Hichem Snoussi

Paul Honeine

DISTRIBUTED LOCALIZATION IN WIRELESS SENSOR NETWORKS AS A PRE-IMAGE PROBLEM IN A REPRODUCING KERNEL HILBERT SPACE

In this paper, we introduce a distributed strategy for localization in a wireless sensor network composed of limited range sensors. The proposed distributed algorithm provides sensor position estimation from local similarity measurements. Incremental Kernel Principal Component Analysis techniques are used to build the nonlinear manifold linking anchor nodes. Non-anchor nodes positions are estimated by the pre-image of their nonlinear projection onto this manifold. This non-linear strategy provides a great accuracy when data of interest are highly corrupted by noise and when sensors are not able to estimate their Euclidean inter-distances.

INTRODUCTION

Recent technological advances in electronics and wireless communications have led to the development of tiny, low-power and low-cost sensors for physical observations purpose. Deployed randomly and densely in the environment of interest, and designed with efficient distributed algorithm, sensor networks seem to offer several opportunities, specially in monitoring and tracking applications [START_REF] Estrin | Instrumenting the world with wireless sensor networks[END_REF]. The first step of estimating the sensorlocation after deployment is thus a crucial issue. GPS system may solve in practice localization problem for each node of the embedded network. However, GPS receivers at each device may be too expensive and too power-intensive for the desired application, whose low energy consumption is the main constraint to respect. As a consequence, we consider only few sensors, called anchor nodes, which have a perfect a priori knowledge of their coordinates thanks to GPS receivers.

The most famous techniques used for localization applications are based upon either semidefinite programming (SDP) or multidimensional scaling (MDS) algorithms (see [START_REF] Bachrach | Localization in sensor networks[END_REF] and references therein). In the SDP framework, proximity measurements between sensors are expressed as geometric constraints, leading to a convex optimization problem. Unfortunately, this approach is unable to accommodate precise range data and unadapted to large scale sensor networks. The MDS algorithm [START_REF] Cox | Multidimensional Scaling, Second Edition[END_REF], which is strongly related to the linear Principal component Analysis (PCA), has shown its effectiveness to estimate unknown sensors location by applying an orthogonal basis transformation. However, if the data are not inter-sensor distances or are linked to coordinates by an unknown non-linear function, linear techniques such as MDS and PCA fail to accurately estimate the node positions.

In this paper, we propose a distributed strategy for solving the localization problem, by borrowing stateof-art methods from machine learning. Each step of the proposed distributed scheme is built in the context of wireless sensor networks application, i.e keeping in mind energy reserve and computational limitations. We assume that each device determines non-Euclidean similarity measurements with other sensors, from some measurements such as the received signal strength indication (RSSI) or estimated covariance sensor data [START_REF] Patwari | Manifold learning algorithms for localization in wireless sensor networks[END_REF]. Thus, we propose to use these similarities [δ ij ] N i,j=1 between neighbor nodes for location estimation purpose. The main idea is to design a nonlinear manifold via a high-dimensional Reproducing Kernel Hilbert Space H (RKHS) thanks to similarities between anchor-node measurements. Next, non-anchor sensors coordinates are estimated by the pre-image of their projection onto the manifold. This paper is organized as follows: Section 2 is devoted to the main contribution of this work which consists of a kernel-based non-linear method for selflocalization. In section 3, simulation results confirming the algorithm efficiency are shown.

Problem statement

Consider a network of N sensor nodes, of m anchors of known positions and N -m unaware-position sensors, living in a p-dimensional space (p = 2 for localization in plane, with N -m m > p). Let x i ∈ IR p be the i-th sensor coordinates, {x i } m i=1 is the set of anchors nodes coordinates whose positions are known, {x i } N i=m+1 fit the unknown remaining sensor coordinates. If we assume that maximum spotting sensor range is equal to a distance r, sensor i will consider sensor j as a neighbor when the distance x ix j is lower than r. The i-th sensor neighborhood is denoted by V(i).

DISTRIBUTED LOCALIZATION ALGORITHM

The proposed algorithm is implemented in three steps. The first step is devoted to build the Hilbert space H associated to the reproducing positive definite kernel κ(x i , x j ) which best approximates estimated anchor pairwise similarities δ ij . The second step is dedicated to map the data into a high-dimensional feature space, obtained from anchor informations by a Kernel-PCA technique. The third part is aimed to reconstruct the N -m unknown sensors positions by a pre-image optimization

Inputs: {[δ ij ] j∈V(i) } m i=1 , {x i } m i=1 Initialization: σ * ← 0 for i = 1 to m • Compute σ * i maximum of Ω i (x 1 , . . . , x m , K * ) at anchor number i • σ * ← σ * + 1 m σ * i
• Communicate σ * to anchor number i + 1 end for Communicate σ * to the N sensors Table 1: Pseudo-code of distributed alignment maximization implemented on the m anchor-nodes, where each anchor knows only nearby anchor positions scheme. The choice of Kernel-PCA [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF] method is supported by its non-linearity property, its flexibility with a large variety of kernel functions and its distributed capabilities.

Kernel selection from anchor similarities

In order to build a valid RKHS, a kernel function respecting the anchor similarities should be selected thanks to the alignment method. The alignment criterion is a measure of similarity between two reproducing kernels or between a kernel and a target function [START_REF] Cristianini | On kernel target alignment[END_REF]. Let K * be a target matrix, and K σ the Gram matrix associated to the reproducing kernel κ σ , with a tuning parameter σ, for a given training set, i.e. with entries κ σ (x i , x j ) for i, j ∈ {1, . . . , m}. The alignment between these two matrices is defined by

A(K σ , K * ) = K σ , K * F K σ , K σ F K * , K * F (1)
where . , . F is the Frobenius inner product between two matrices. For our purpose, the target matrix K * is given by our similarity matrix, with

K * (i, j) = δ ij , (2)
for all i, j ∈ {1, . . . m}, in which case K * , K * F is constant. We have restricted potential kernels to Gaussian kernels defined as:

κ σ (x i , x j ) = exp(-x i -x j 2 /2σ
2 ) where σ is a positive scalar to be determined and x i , x j are anchor coordinates. Therefore, the aim is to maximize the alignment criterion (1) with respect to σ.

Thanks to Lagrange method, the optimization problem is equivalent to:

σ * = arg max σ m i=1 m j=1 (δ ij κ σ (x i , x j )) -λ m i,j=1 κ σ (x i , x j ) 2 ,
where the term between brackets will be denoted by

Ω i (x 1 , . . . , x m , K * ).
In practice, if the anchor j is a neighbor of anchor i, δ ij is computed, otherwise it is set to zero. This formulation respects the hypothesis of spatial correlation [START_REF] Patwari | Manifold learning algorithms for localization in wireless sensor networks[END_REF] and helps to preserve energy by limiting communications to neighbors. Thus, we obtain

Ω i (x 1 , . . . , x m , K * ) = j∈V(i) δ ij κ σ (x i , x j )-λ j∈V(i) κ σ (x i , x j ) 2
Pre-processing: where λ is the Lagrange coefficient and V(i) is the set of anchor-neighbors of anchor i. A distributed gradient algorithm descent is performed from anchor to anchor according to the scheme illustrated in Table 1.

learn σ * as in Table 1 Inputs: K = [κ σ * (x i , x j )] m i,j=1 Initialization: Randomly set A(0) for k = 1 to m • Compute A(k) from (4)
After computing σ * as the best alignment, the second step consists in building the nonlinear manifold related to a subspace in the RKHS H.

Kernel-PCA upon anchor nodes

As a nonlinear extension of PCA, Kernel-PCA determines the principal axes in a RKHS, which can be constructed explicitly from the input space by a nonlinear map φ(.), or implicitly by considering the corresponding reproducing kernel κ(•, •). In other words, for an m input data problem, [x k ] m k=1 , a PCA is performed in H yielding a set of m -1 orthogonal axes 1 {v 1 , . . . , v m-1 }, and thus defining a subspace P of H. Since the latter lies in the span of the φ-images of the input data, Kernel-PCA [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF] is computed by diagonalizing the dot-product matrix K in H, by solving

mλ k a k = Ka k 1 ≤ k ≤ m -1 (3)
for the column vectors a k , with k = 1, . . . , m -1, of feature coefficients, and K = [ φ(x i ), φ(x j ) ] m i,j=1 is the dot-product matrix in H, known as the Gram matrix or learning matrix. This is done without the need to carry out the map φ explicitly. Only the kernel values κ(x i , x j ) are needed.

For the choice of the reproducing kernel, there is no restrictions with Kernel-PCA. We propose to use the one obtained from section 2.1, with the maximum alignment criterion. Therefore, we consider the Gram matrix [κ σ * (x i , x j )] m i,j=1 . The mapping φ : IR p → H corresponds to a RKHS H where dot products φ(x i ), φ(x j ) H are as close as possible to similarity measurements δ ij .

We assume that each anchor-sensor has at least one anchor as a neighbor node which means that all anchor nodes respect a connected graph. As considered in section 2.1 and illustrated in Table 1, each anchor already knows the positions of its neighbor anchors, which for sensor i is given by {x j } j∈V(i) . Thus, it can compute κ σ * (x i , x j ) for j ∈ V(i) or nullify it otherwise. 1 For the sake of clarity, we assume that all mapped training pattern are centered in the feature space. Therefore, we count m -1 eigenvectors associated to positive eigenvalues non-null.

Otherwise one should substitute K in (3) with K -1 m 1mK -

1 m K1m + 1 m 2 1mK1m where (1m) ij = 1 placements x a x b x c x i xi φ(.) φ(x i ) φ(x a ) φ(x b ) φ(x c ) P P φ(x i )
Figure 1: Illustration of the idea behind Kernel-PCA and pre-image techniques for localization in sensor networks. The left-hand-side frame corresponds to the input space (plane) and the right-hand-side to the feature space. The latter is reduced to the subspace P defined by the anchor nodes (red squares). The position of any non-anchor node (black circle) is estimated from ( 7) by a pre-image technique, after projecting its image onto P.

After evaluating the anchor-node pairwise similarities, an incremental Kernel-PCA algorithm is run over the m anchors of the network by considering the kernel values [κ σ * (x i , x j )] m i,j=1 . The Kernel Hebbian algorithm [START_REF] Kim | Iterative kernel principal component analysis for image modeling[END_REF], presented as a direct application of the Generalized Hebbian algorithm in an RKHS, is dedicated to solving a kernel eigen-problem in a distributed way. The matrix of feature coefficients

A = [a 1 a 2 • • • a m-1 ] ∈ IR m×m-1
is updated by each anchor node as follows:

A(k + 1) = A(k) + η y(k)b(k) t -LT y(k)y(k) t A(k) , (4)
where η is a predefined learning rate,

y(k) = m j=1 a j (k)κ σ * (x k , x j ), b(k) = [0 • • • 1 • • • 0]
a canonical unit vector whose k -th element is 1, and LT (•) the lower triangular operator.

At the initialization stage, A(0) is assigned a random value2 by the anchor number 1. Then it is communicated between anchors and updated with respect to the recursive expression (4). After this learning stage, A * is communicated to the remaining N -m sensors of the network. The pseudo-code is described in Table 2. Next, for each sensor i, we represent its image φ(x i ) in the space defined by the anchors, and obtained from the feature vectors defined by A * . The problem of estimating xi from that representation is the classical pre-image reconstruction problem.

Pre-image for location estimation

As we consider a reproducing kernel associated to a nonlinear map φ(•), the induced RKHS is of higher dimensionality, and infinite for some kernels such as the Gaussian kernel. Thus, the image φ(x i ) of any sensor i has redundant information. We consider representing it with the finite-number features in A * obtained with anchors, as we hope that the remaining dimensions contain only noise information. Representing the data in an eigen-space is known in machine learning literature as the empirical kernel map [START_REF] Schölkopf | Learning with Kernels[END_REF]. Next we study the problem of getting back to the input space, with a pre-image technique, in order to estimate the N -m coordinates of non-anchor sensors. But before, we recall that P is spanned by the m -1 eigenvectors, or less if some eigenvalues are null, and define P as the projection operator onto it. Thus, P φ(x) is the projection of φ(x) onto P, as illustrated in Fig. 1.

The projection satisfies the minimal reconstruction error P φ(x) -φ(x) 2 . While we have no access to x, we want to determine x ∈ IR p whose image φ(x) best approximates P φ(x), in the least square sense. Since P φ(x) is already known given only the anchor/nonanchor similarities, the optimized problem could then be put in the following form [START_REF] Mika | Kernel PCA and de-noising in feature spaces[END_REF]:

x = arg min z P φ(x) -φ(z) 2 (5)
This optimization problem should be solved at each of the N -m non-anchor sensors, in order to estimate their unknown coordinates. Each non-anchor sensor j has to evaluate m kernel values κ σ * (x i , x j ) for i ∈ {1, . . . , m} (i.e. with the m anchors). However, this computational cost is reduced in practice, as a small number of anchors are in the neighborhood of each sensor, i.e if anchor i is inside the visibility range of sensor j, we have

κ σ * (x i , x j ) = δ ij , otherwise κ σ * (x i , x j ) = 0.
We remind that non-anchor sensors have already received the last update of expansion matrix A * . Each sensor can extract its nonlinear principal components. Assuming that β k j is the projection onto the k-th component for the non-anchor sensor j, it reads:

β k j = m i=1 a i k κ σ * (x i , x j ) (6)
Therefore, we can write P φ(x j ) = m-1 k=1 β k j v k where x j is the sensor j unknown coordinates, β k j is the pro-Pre-processing:

learn A * as in Table 2 Inputs: A * , σ * , {[δ ij ] l∈V(j) } N j=m+1 , {x i } m i=1
Initialization: randomly set x m+1 , . . . , x N for j = m + 1 to N • At node number j, compute xj from (7) end for jection on the k th component and v k is the k-th eigenvector. Solution from ( 5) can be simplified in terms of dot products in H:

xj = arg min z κ σ * (z, z) -2P φ(x j ) t .φ(z)
where we neglected a term independent of z. If the considered kernel is radial such as the Gaussian kernel, thus κ(z, z) is constant for all z, we get the following optimization problem (see [START_REF] Mika | Kernel PCA and de-noising in feature spaces[END_REF] for more details)

xj = arg max z P φ(x j ) t φ(z) = arg max z m-1 k=1 β k j m i=1 a i k κ σ * (z, x i ), (7)
Thus, a gradient descent algorithm is performed by nonanchor sensor j in order to minimize the criterion [START_REF] Kim | Iterative kernel principal component analysis for image modeling[END_REF]. Initially, all non-anchors nodes have a random position in the network and the optimization is executed locally with no synchronization constraints. The preimage technique for localization in sensor networks is illustrated in Table 3.

EXPERIMENTS

The algorithm proposed in this paper is independent of the type of similarities considered for localization purpose. In our simulation scenario, we make the assumption that data inputs are jointly distributed with a covariance being function of the distance. Thus, estimated covariance relations are used as local similarity measurements between neighbor nodes [START_REF] Patwari | Manifold learning algorithms for localization in wireless sensor networks[END_REF]. For this, let a network consists of sensors measuring the same physical phenomena, such as temperature, pressure or luminance measurements. Here, w i represents the vector data recorded by sensor i on a given time interval. By considering a static field, we assume that data readings {w i } N i=1 are jointly generated from a normal distribution of mean µ = [µ 1 • • • µ N ] and covariance matrix C shaped as where ψ : [0, ∞[→ [0, 1] is a non-negative decreasing function. The experimental setting consists of 20 anchor-sensors and 80 non-anchor sensors randomly spread over a square surface 100 × 100. The data readings for each sensor consists of 200 measurements. We consider the Gaussian kernel, while its bandwidth σ * is given by maximizing the alignment.

C = [ψ( x i -x j )] N i,j=1 (8)
In the first experience (noted Exp.I), data are generated according to the unknown covariance ψ 1 (z) = exp(z 2 2τ 2 ) where τ = 20. Localization results are shown in Fig. 5(a) with sensor range r = 45. Fifty simulations were run for each range value. Results are presented in Fig. 2 for different sensor range value, the Lagrange coefficient for alignment λ is fixed to 0.3.

The second experiment (noted Exp.II) is dedicated to test the robustness of our algorithm by using a generative covariance ψ different from the Gaussian kernel. The spherical model [START_REF] Gneiting | Compactly supported correlation functions[END_REF], which is commonly applied in environmental and geological sciences, is used with where d is a cut-off distance. In Exp.II, the parameter d is fixed to 60 in order to obtain a decreasing speed of covariance with distance close to that of Exp.I with τ = 20 (see Fig. 4). An example of localization results is shown in Fig. 5(b) with r = 45. Note that some similarity measurements are missed when the sensor range r is fixed to 45 and τ = 20 or d = 60 (see Fig. 5), since remote sensor pairwise are ignored in each algorithm step.

ψ 2 (z) = 1 -3 2d z + 1 2d 3 z 3 for 0 ≤ z ≤ d 0 for d < z
Results of fifty simulations with λ = 0.3, are shown in Fig. 3 for different sensor range values.

In both experiments, we note that localization performances are significantly degraded with low range (around a range of 20). However, localization error quickly decreases when number of neighbor anchors rises above 3. This is mainly due to a constraint in local localization, since the subproblems are well defined when 3 anchors available in IR 2 . We note that localization error is always sensitively better in Exp.I than in Exp.II. This is justified by the similarity between the considered Gaussian kernel and the generative covariance function ψ 1 . Nevertheless, results of Exp.II show that our algorithm is still robust when faced to more realistic situations where no similarity exists between the chosen Gaussian kernel and the underlying generative covariance.

CONCLUSION AND PERSPECTIVES

In this article, we have shown that Kernel-PCA, coupled with reconstruction techniques, computes unknown coordinates of the network sensors thanks to a non-linear transformation φ of the input space. This localization problem is solved distributively with limited sensor range and noisy non-Euclidean measurements. Good performances in positioning are observed thanks to the non-linear process. There are several directions for further research including a way to add additional information in pre-image criterion, such as similarity measurements between neighbor non-anchor nodes.

Figure 2 :

 2 Figure 2: Evolution of mean distance error per sensor with a gaussian generative covariance function τ = 20 (Exp.I).

Figure 3 :

 3 Figure 3: Evolution of mean distance error per sensor with a third degree polynomial generative covariance function d = 60 (Exp.II).

Figure 4 :

 4 Figure 4: Graph comparing curves of functions ψ1 with parameter τ = 20 and ψ2 with d = 60.

Figure 5 :

 5 Figure 5: The 20 anchors are represented by blue squares. Estimated location of 80 nodes are indicated by a red cross, distance error in location to true position are symbolized by segment. (a) and (b) are respectively related to Exp.I and Exp.II for the same initial configuration of sensors and range fixed to 45.

Table 2 :

 2 at anchor k • Communicate A(k) to anchor k + 1 end for Communicate A * to the non-anchor nodes Pseudo-code of incremental Kernel-PCA algorithm implemented on the m anchor-nodes.

Table 3 :

 3 Pseudo-code of pre-image algorithm for localization estimation.

A(0) must neither be the zero vector nor orthogonal to the eigenvectors.