Jie Chen
email: jie.chen@utt.fr

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

Paul Honeine
email: paul.honeine@utt.fr

Hichem Snoussi
email: hichem.snoussi@utt.fr

Henri Lantéri
email: henri.lanteri@unice.fr

Céline Theys
email: celine.theys@unice.fr

Techniques d'apprentissage non-linéaires en ligne avec contraintes de positivité

Keywords: apprentissage, régression, non-linéaire, positivité, adaptatif, réseau de capteurs

Cet article décrit une nouvelle classe d'algorithmes d'apprentissage non-linéaires en ligne avec contrainte de positivité sur la solution. Ceux-ci sont appliqués au problème d'identification distribuée d'un champ scalaire positif, par exemple de rayonnement thermique ou de concentration d'une espèce chimique, par un réseau de capteurs. La question du suivi de l'évolution de la grandeur physique surveillée au cours du temps est également considérée. Les algorithmes proposés sont testés sur des données synthétiques régies par des équations de diffusion. Ils démontrent une excellente capacité de suivi des évolutions du système, tout en affichant un coût calculatoire réduit.

I. Introduction

Le domaine des réseaux de capteurs sans fil fait actuellement l'objet d'un intérêt considérable de la part des communautés académique et industrielle. La dispersion d'une multitude de capteurs bon marché dans une région donnée, l'élaboration d'un protocole de routage adéquat, et une implémentation algorithmique efficace, ouvrent en effet de nombreuses perspectives d'applications civiles et militaires. Dans cet article, on s'intéresse au problème d'identification distribuée d'un champ scalaire positif, par exemple un rayonnement thermique ou la concentration d'une espèce chimique, et le suivi de son évolution au cours du temps.

Le mode de calcul distribué est inhérent au caractère réparti des noeuds du réseau sur l'aire surveillée, dont la tâche est d'acquérir et traiter localement les mesures. L'efficacité de la procédure d'identification est conditionnée par les interactions des noeuds, dictées par la topologie du réseau. Deux principes de coopération ont été principalement considérés dans la littérature. En mode incrémental, l'information transite de façon séquentielle et cyclique d'un noeud voisin à l'autre. Le coût énergétique, largement dicté par le volume des communications, tend à être minimum [START_REF] Lopes | Incremental adaptive strategies over distributed networks[END_REF], [START_REF] Rabbat | Distributed optimization in sensor networks[END_REF]. En mode de diffusion, chaque noeud coopère avec l'ensemble de ses voisins pour une meilleure qualité d'estimation [START_REF] Lopes | Diffusion least-mean squares over adaptive networks[END_REF]. Il en résulte un nombre accru d'échanges d'informations, qu'il est possible de restreindre en limitant artificiellement le nombre de voisins concernés. On parle dans ce cas de diffusion probabiliste.

La complexité des applications considérées ici, telles que la diffusion de chaleur mesurée par rayonnement thermique qui illustre la suite de cet article, nécessite de recourir à des méthodes d'identification adéquates. Fondées sur les travaux précurseurs d'Aronszajn [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], les récentes avancées de l'estimation fonctionnelle basée sur les espaces de Hilbert à noyau reproduisant, et de la théorie de la régularisation, ont apporté des réponses convaincantes aux problèmes implicant des systèmes non-linéaires. La stratégie d'apprentissage en ligne considérée, reposant sur un mode coopératif de type incrémental, démontre d'excellentes capacités de suivi des évolutions du système tout en affichant un coût calculatoire réduit. Si plusieurs approches similaires existent dans la littérature, l'article présenté tire en partie son originalité de l'usage du formalisme des espaces de Hilbert à noyau reproduisant.

Ce travail fait suite à [START_REF] Honeine | Distributed regression in sensor networks with a reduced-order kernel model[END_REF], [START_REF] Honeine | Functional estimation in Hilbert space for distributed learning in wireless sensor networks[END_REF] dont il se distingue par la contrainte de positivité pesant sur le champ scalaire estimé, ce qui constitue une originalité par rapport à l'ensemble de la littérature du domaine. On s'inspire pour cela de méthodes de déconvolution, où la prise en compte d'informations a priori sur l'espace des solutions admissibles est primordiale pour garantir la stabilité du résultat et traduire le phénomène physique associé. Les algorithmes multiplicatifs jouent un rôle central dans la résolution de problèmes inverses sous contrainte de positivité, en raison de leur formulation simple assurant très naturellement la non-négativité de la solution. Les plus populaires sont l'algorithme de Richardson-Lucy [START_REF] Lucy | An iterative technique for the rectification of observed distributions[END_REF], [START_REF] Richardson | Bayesian based iterative method of image restauration[END_REF], obtenu par maximisation de la vraisemblance sous hypothèse poissonienne, particulièrement utilisé en imagerie astrophysique et médicale, et l'algorithme ISRA de Daude-Witterspoon et Muehlenner [START_REF] Daude-Witherspoon | An iterative image space reconstruction algorithm suitable for volume ECT[END_REF], adapté à des données gaussiennes. La méthode SGM -Split Gradient Method -de Lantéri et coll. [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF] généralise ces techniques en les groupant sous une formulation unique, dédiée à la minimisation de divergences entre données et modèles sous contrainte de positivité. L'article est organisé ainsi. On adapte d'abord l'algorithme SGM au problème d'identification considéré. Son implémentation est distribuée de sorte à satisfaire aux contraintes d'un réseau de capteurs coopérant en mode incrémental, et ses capacités de poursuite développées par un calcul de gradient stochastique. On rappelle ensuite les fondements de l'estimation fonctionnelle dans les espaces de Hilbert à noyau reproduisant, que l'on traduit algorithmiquement afin d'être à même de traiter des phénomènes non-linéaires complexes. L'approche obtenue est enfin mise en oeuvre sur des données synthétiques. Une conclusion et des perspectives viennent clore la discussion.

II. Formulation générale de la méthode

Etant donné un ensemble de couples (x i , y i ) ∈ IR p + ×IR + , on s'intéresse au problème d'identification d'un modèle de la forme générale y = ψ(x), où ψ appartient à un espace fonctionnel H défini ultérieurement. Afin de fixer les idées, on considère dans un premier temps un modèle linéaire y i = α x i . Pour de nombreuses applications, une contrainte de non-négativité sur les composantes de α s'impose naturellement. Voir [START_REF] Fevotte | Nonnegative matrix factorization with the Itakura-Saito divergence. Application to Music Analysis[END_REF], [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF]. Celle-ci est au coeur de la méthode développée ici, appliquée en fin d'article à l'estimation d'un champ scalaire positif par un réseau de capteurs. Le problème à résoudre se formule ainsi

α o = arg min α J(α) (1)
sous contrainte de α ≥ 0 (2) où J(α) désigne une fonction coût donnée. Le lagrangien correspondant s'écrit

Q(α, λ) = J(α) -λ α (3)
avec λ le vecteur constitué des facteurs de Lagrange, tous non-négatifs. Les conditions de Kuhn-Tucker à l'optimum se traduisent par

∇ α Q(α o , λ o ) = 0 (4) λ o i α o i = 0 ∀i, (5)
avec α o i (resp. λ o i ) la i-ème composante de α o (resp. λ o ). Ceci se résume par le problème suivant à résoudre

α o i [∇J(α)] i = 0, (6)
où [∇J(α)] i désigne la i-ème composante de ∇J(α). Pour aboutir à un algorithme d'optimisation de type multiplicatif décrit dans la suite, les auteurs de [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF] suggèrent d'appliquer la méthode du point fixe à [START_REF] Honeine | Distributed regression in sensor networks with a reduced-order kernel model[END_REF]. Celle-ci mène à 1 α

(k+1) i = α (k) i + η (k) i f i (α (k) ) α (k) i [-∇J(α (k) )] i , (7) avec η
(k) i > 0 un facteur de relaxation et f i (α) une fonction positive sur le domaine admissible de α. Ces derniers et le signe négatif introduit dans l'expression ci-dessus sont respectivement destinés, d'une part à pouvoir agir localement sur le caractère contractant de la fonction ainsi traitée par la méthode du point fixe, d'autre part à adopter pour direction de descente celle qui est opposée au gradient.

Garantir la positivité de α

(k+1) i à partir de celle de α (k) i impose la condition η (k) i ≤ 1 f i (α (k) ) [∇J(α)] i (8)
si [∇J(α)] i > 0. Dans le cas contraire, où [∇J(α)] i ≤ 0, aucune restriction relative à la contrainte de positivité ne s'impose dans le choix de ce pas. On en déduit finalement l'expression générale de l'algorithme

α (k+1) = α (k) + η (k) d (k) , (9)
avec d (k) la direction de descente définie par

d (k) = -diag[f i (α (k) ) α (k) i ] ∇J(α (k) ), (10)
et η (k) le pas résultant éventuellement d'un algorithme de recherche linéaire [START_REF] Bertsekas | Nonlinear Programming[END_REF] dans l'intervalle ]0, η

(k) max ] avec η (k) max = min i 1 f i (α (k) ) [∇J(α (k) )] i . ( 11
)
On s'intéresse à présent à l'algorithme multiplicatif associé [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF]. Cette forme est récemment devenue très populaire en raison notamment de sa simplicité de mise en oeuvre. Elle sera exploitée par la suite dans le cadre des réseaux de capteurs. Pour ce faire, on décompose l'opposé du gradient -∇J(α (k) ) ainsi

-∇J(α (k) ) = U (α (k) ) -V (α (k) ), (12)
où U et V sont deux fonctions à valeurs strictement positives. Voir [START_REF] Lantéri | Méthode algorithmique de minimisation de fonctions d'écart entre champs de données. Application à la reconstruction d'images astrophysiques[END_REF] pour une discussion sur la non-unicité de cette décomposition et ses interprétations. En reportant celle-ci dans l'expression [START_REF] Honeine | Functional estimation in Hilbert space for distributed learning in wireless sensor networks[END_REF], on obtient

α (k+1) i = α (k) i +η (k) i f i (α (k) ) α (k) i [U (α (k) )-V (α (k) )] i . (13)
En posant

f i (α (k) ) = 1/[V (α (k) )] i , on aboutit à α (k+1) i = α (k) i + η (k) i α (k) i [U (α (k) )/V (α (k) ) -1] i . ( 14
)
Comme dans [START_REF] Lopes | Diffusion least-mean squares over adaptive networks[END_REF], garantir la positivité de α

(k+1) i impose la condition η (k) i ≤ 1 1 -[U (α (k) ] i /[V (α (k) )] i (15) si [∇J(α)] i > 0, c'est-à-dire [V (α (k) )] i > [U (α (k) )] i .
Dans le cas contraire, aucune condition ne s'impose. On note que le second membre de l'inégalité [START_REF] Richardson | Bayesian based iterative method of image restauration[END_REF] est strictement supérieur à 1. En prenant un pas unitaire dans [START_REF] Richard | Online prediction of time series data with kernels[END_REF], qui garantit en conséquence la positivité du résultat, on aboutit à la formulation multiplicative suivante

α (k+1) i = α (k) i [U (α (k) )] i /[V (α (k) )] i . (16)
En notant ⊗ (resp. ) le produit (resp. la division) au sens d'Hadamard, on reformule enfin ce résultat ainsi

α (k+1) = α (k) ⊗ U (α (k) ) V (α (k) ). ( 17
)
Il convient de noter le peu d'opérations mises en jeu par cet algorithme, qui assure naturellement la positivité de la solution tout au long du processus d'optimisation si la condition initiale l'est. On relève également que les composantes nulles de α (k) se propagent au cours des itérations, ce qui peut constituer un avantage dans la recherche de solutions parcimonieuses. Enfin, on constate que la formulation multiplicative (17) ne repose sur aucune expression particulière du critère J(α) utilisé, ni du modèle sous-jacent.

III. Description de l'algorithme distribué Dans le contexte de l'apprentissage distribué dans un réseau de N capteurs, on s'intéresse au problème d'identification d'un champ scalaire positif, par exemple un rayonnement thermique ou la concentration d'une espèce chimique, et le suivi de son évolution au cours du temps. Chaque noeud n acquière pour cela des mesures (x n, , y n, ) au cours du temps. Pour simplifier l'exposé, on suppose que le critère J(α) peut être décomposé en une somme de critères individuels J n (α), un en chaque noeud n, soit

J(α) = N n=1 J n (α). ( 18
)
Pour déterminer la solution α o , une approche naturelle consiste à mettre en oeuvre l'algorithme de descente ( 9) sous la forme

α (k+1) = α (k) + N n=1 η (k) d (k) n , (19)
avec d (k) n la direction de descente donnée ici par

d (k) n = -diag f i (α (k) ) α (k) i ∇J n (α (k) ), ( 20
) et η (k) un pas choisi dans l'intervalle ]0, η (k)
max ] avec

η (k) max = min i 1 f i (α (k) ) N n=1 ∇J n (α (k) ) i . (21)
Il en résulte l'algorithme incrémental provisoire suivant où les paramètres du modèle sont mis à jour, successivement, d'un noeud voisin à l'autre. 1.

β (k) 0 = α (k) 2. β (k) n = β (k) n-1 + η (k) d (k) n , n = 1, . . . , N 3. α (k+1) = β (k) N
Cette approche nécessite cependant, par l'instruction 2., que chaque noeud ait accès à l'ensemble des données et à α (k) pour le calcul de η (k) et d (k) n . Cet inconvénient se retrouve dans la formulation multiplicative qui, après qu'on ait posé -∇J n (α

(k) ) = U n (α (k) ) -V n (α (k) ), s'écrit donc α (k+1) = α (k) ⊗ N n=1 U n (α (k) ) N n=1 V n (α (k) ). ( 22
)
Afin de résoudre ce problème, qui entraînerait un coût de communication prohibitif, on préconise une évaluation locale du pas et de la direction de descente au point reçu du précédent capteur. L'équation de mise à jour est alors donnée par

β (k) n = β (k) n-1 + η (k) n d (k) n , (23)
avec d (k) n la direction de descente définie par

d (k) n = -diag f i (β (k) n-1 ) β (k) n-1,i ∇J n (β (k) n-1 ), (24) et η (k)
n un pas choisi dans l'intervalle ]0, η

max ] avec

η (k) max = min i 1 f i (β (k) n-1 ) ∇J n (β (k) n-1 ) i . (25)
L'algorithme itératif proposé, élaboré sur le mode de coopération incrémental, est finalement donné par Algorithme itératif

A chaque instant k, répéter 1. β (k) 0 = α (k) 2. β (k) n = β (k) n-1 + η (k) n d (k) n , n = 1, . . . , N 3. α (k+1) = β (k) N
Le réseau est parcouru séquentiellement, et les paramètres du modèle transmis d'un noeud voisin à l'autre, afin de suivre l'évolution du phénomène physique au cours temps.

En posant f i (β

(k) n-1 ) = 1/[V n (β (k) n-1 )] i et η (k) n
= 1, on aboutit finalement à l'algorithme multiplicatif suivant.

Algorithme multiplicatif

A chaque instant k, répéter 1. β (k) 0 = α (k) 2. β (k) n = β (k) n-1 ⊗ U n (β (k) n-1 ) V n (β (k) n-1 ), n = 1, . . . , N 3. α (k+1) = β (k) N
La section suivante vise à préciser ces algorithmes sous des hypothèses statistiques particulières, par une sélection appropriée de la fonction coût utilisée.

IV. Choix de fonctions coût

Il convient de préciser maintenant quelques choix possibles pour la fonction coût J(α). On motive ceux-ci par des hypothèses statistiques, en supposant dans un premier temps que les données sont poissoniennes, en considérant dans un second temps qu'elles sont entachées d'un bruit gaussien additif. Pour autant, ces choix ne se limitent pas à ceux développés ci-dessous.

A. Cas de données poissoniennes

Dans le cadre d'un procédé parfait de détection de photons par un réseau de capteurs, destiné à la mesure d'un rayonnement thermique par exemple, l'intensité y n, en chaque capteur n est supposée suivre une loi de Poisson. On fait l'hypothèse que la moyenne de celle-ci est régie par un modèle linéaire y n, = α x n, , où (x n, , y n, ) est un couple de données acquises par le capteur n, supposées non-négatives. A titre de généralisation, on suppose que L n de ces couples de données sont disponibles au moment de la mise à jour des paramètre du modèle par le capteur n. La fonction de vraisemblance s'écrit aisément

F (α) = N n=1 Ln =1 α x n, y n, y n, ! exp(-α x n, ) (26)
Par la formule de Stirling, la fonction de log-vraisemblance s'exprime selon

-log F (α) ≈ N n=1 Ln =1 (α x n, -y n, ) + y n, log y n, α x n, . (27)
Il s'agit là d'une I-divergence de Csiszär entre α x n, et y n, . Cette mesure généralise la divergence de Kullback à des fonctions non-négatives qui ne sont pas des distributions. Le problème à résoudre se formule ainsi

α o = arg min α N n=1 J n (α) (28) sous contrainte que α ≥ 0 (29) où J n (α) = Ln =1 (α x n,
) -y n, log(α x n, ) désigne la fonction coût associée au capteur n. Les résultats qui suivent découlent directement des développements de la Section III. En l'occurrence, il convient de calculer le gradient de J n (α) au préalable

∇J n (α) = Ln =1 x n, - y n, α x n, x n, . (30)
Si la décomposition de -∇J n (α) n'est pas unique, il vient immédiatement que l'on peut toutefois poser

U n (α) = Ln =1 y n, α x n, x n, V n (α) = Ln =1 x n, . (31)
L'algorithme correspondant, sous sa forme multiplicative, s'écrit finalement ainsi Algorithme multiplicatif (I-divergence de Csiszär)

A chaque instant k, répéter 1. β (k) 0 = α (k) 2. β (k) n = β (k) n-1 ⊗ Ln =1 y n, x n, β (k) n-1 x n, Ln =1 x n, 3. α (k+1) = β (k) N
Il convient de noter la simplicité de cet algorithme, qui vise pourtant à minimiser une I-divergence de Csiszär en assurant la positivité de la solution tout au long du processus d'optimisation.

B. Cas de données gaussiennes

On suppose que l'intensité y n,k en chaque capteur n est ici encore régie par un modèle linéaire y n, = α x n, , à présent corrompue par un bruit additif gaussien i.i.d. centré de variance σ 2 n . La fonction de vraisemblance est

F (α) = N n=1 Ln =1 1 2πσ 2 n exp - (α x n, -y n, ) 2 2σ 2 n . (32)
La fonction de log-vraisemblance correspondante s'écrit dans ce cas

-log F (α) = 1 2 N n=1 Ln =1 (α x n, -y n, ) 2 σ 2 n , (33)
d'où l'on a écarter les termes constants additifs. Comme précédemment, le problème à résoudre se formule ainsi

α o = arg min α N n=1 J n (α) (34) sous contrainte que α ≥ 0 (35) où J n (α) = 1 2σ 2 n k (α x n,k -y n,k
) 2 désigne la fonction coût associée au capteur n. Le calcul du gradient de cette dernière conduit à

∇J n (α) = 1 σ 2 n Ln =1 (x n, x n, α -y n, x n, ). (36)
Il apparaît immédiatement que la décomposition suivante de -∇J n (α) en une différence de fonctions strictement positives convient

U n (α) = 1 σ 2 n Ln =1 y n, x n, ( 37
) V n (α) = 1 σ 2 n Ln =1 x n, x n, α. (38)
On en déduit l'algorithme multiplicatif suivant Algorithme multiplicatif (erreur quadratique)

A chaque instant k, répéter 1. β (k) 0 = α (k) 2. β (k) n = β (k) n-1 ⊗ Ln =1 y n, x n, Ln =1 x n, x n, β (k) n-1 3. α (k+1) = β (k) N
Il est intéressant et rassurant de noter que l'algorithme ci-dessus fait intervenir la matrice de corrélation empirique des données d'entrée, ainsi que le vecteur d'intercorrélation empirique entre les entrées et sorties. On constate enfin que les deux algorithmes développés au cours de cette section se confondent lorsqu'une unique donnée est disponible par capteur, au moment de la mise à jour. L'instruction 2. s'écrit alors simplement

β (k) n = y n, x n, β (k) n-1 β (k)
n-1 . (39)
Ces éléments seront approfondis dans un prochain article.

V. Méthodes d'identification non-linéaires La complexité des applications envisagées, telles que la diffusion de chaleur mesurée par rayonnement thermique qui illustre la fin de cet article, nécessite de recourir à des méthodes d'identification adéquates. La présente section vise à dépasser le cadre restrictif des modèles linéaires pour considérer des formes y = ψ(x ) plus générales, où ψ appartient à un espace fonctionnel H qu'il convient de définir attentivement afin de rendre les calculs aisés. Fondées sur les travaux précurseurs d'Aronszajn [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], les avancées de l'estimation fonctionnelle basée sur les espaces de Hilbert à noyau reproduisant, et de la théorie de la régularisation, ont apporté des réponses convaincantes aux problèmes implicant des systèmes non-linéaires. La stratégie d'apprentissage considérée ci-après reposent sur cette théorie.

A. Principe général

On considère un espace de Hilbert H à noyau reproduisant constitué de fonctions réelles sur un compact X ⊂ IR p . On note • , • H son produit scalaire. Soit κ : X × X -→ IR le noyau de cet espace, qui vérifie par conséquent la propriété reproduisante suivante : ψ(x ) = ψ, κ(x , •) H pour toute fonction ψ de H et tout x de X . Cette dernière étant vérifiée par le noyau lui-même, on a évidemment

κ(x i , x j ) = κ(x i , •), κ(x j , •) H , (40)
ce qui justifie l'appellation de noyau reproduisant. Le noyau gaussien, donné par κ(x i , x j ) = exp(-x i -x j 2 /2σ 2 0 ) avec σ 0 la largeur de bande, en est un exemple.

Etant donné un ensemble {(x , y )} =1,...,L de données d'apprentissage i.i.d., on considère le problème de minimisation d'un critère J entre les sorties du modèle ψ(x ) et les réponses désirées y

ψ o = arg min ψ∈H J (ψ(x ), y ) =1,...,L . (41)
Le Théorème de Représentation Généralisé [START_REF] Schölkopf | A generalized Representer Theorem[END_REF] établit que la solution ψ o peut être recherchée dans l'espace engendré par les L fonctions noyau κ(x , •), c'est-à-dire

ψ o (x) = L =1 α o κ(x , x), (42)

B. Critères de parcimonie

On note que le modèle linéaire généralisé (42) compte autant de termes qu'il y a de données disponibles, limitant son intérêt pratique. Il s'agit là d'une difficulté rencontrée de manière récurrente avec les méthodes à noyau reproduisant. La vaste littérature sur le sujet préconise généralement l'usage de modèles approchés à taille réduite

ψ (k) (x) = |D| i=1 α (k) i κ(x ωi , x) (43)
où les éléments x ωi composent un dictionnaire D. Ce dernier est adapté au cours des itérations, à partir d'un critère complémentaire de parcimonie. Voir [START_REF] Richard | Online prediction of time series data with kernels[END_REF] et bibliographie incluse. On s'en tient ici à une règle simple, tout en sachant que de nombreux raffinements sont possibles. Lors de l'acquisition de tout nouveau couple de données (x , y ), on s'attache dans un premier temps à la mise à jour des coefficients constituant le modèle étendu

ψ (k) (x) + κ(x , x) (44)
avec = max{(y -ψ (k) (x ))/κ(x , x ); 0}. On utilise l'un des algorithmes présentés au cours des sections précédentes pour cela, itératif ou multiplicatif, distribué ou non selon le contexte. Dans un second temps, afin de limiter la taille du dictionnaire D, on ne retient finalement que les composantes excédant un seuil ν 0 donné.

VI. Expérimentations

Pour illustrer la pertinence des algorithmes proposés, on s'intéresse à l'estimation d'un champ de température régi par l'équation différentielle suivante

∂Θ(x, t) ∂t -c∇ 2 x Θ(x, t) = Q(x, t)
à partir du rayonnement thermique associé. Dans cette expression, Θ(x, t) est la température dépendant de la position x et du temps t, ∇ 2 x Θ(x, t) l'opérateur spatial de Laplace, et Q(x, t) la quantité de chaleur apportée. On considère une surface rectangulaire, de conductivité thermique c = 0.1, sur laquelle N capteurs sont tirés uniformément sur une grille de taille 21 × 21. Deux sources de chaleur de 200 W sont activées successivement, la première des instants t = 1 à t = 100 située dans le quart-plan inférieur-droit, la seconde de t = 100 à t = 200 située dans le quart-plan supérieur-gauche.

L'objectif de la simulation proposée est, étant donné des mesures d n, telles que

d n, ∼ P(Θ(x n , t )) (45)
avec P(λ) un bruit de Poisson de paramètre λ traduisant un processus de détection de photons par un capteur, d'estimer Θ(x, t) via ψ (t) (x). Les simulations ont été réalisées avec l'algorithme multiplicatif distribué, développé pour les données poissoniennes, associé à un noyau gaussien de largeur de bande σ égal à 0.18. Le seuil ν 0 destiné à sélectionner les éléments les plus significatifs pour le dictionnaire a été fixé à max i (α i )/10.

Comme l'indique la figure 1, l'algorithme incrémental suit parfaitement les évolutions de la distribution de chaleur. Le choix des capteurs représentés dans le modèle via le dictionnaire D, cerclés de rouge sur chaque représentation, est approprié. On rappelle que leur nombre, ici restreint, correspond à la taille du modèle. La convergence de l'erreur quadratique moyenne normalisée, définie par est illustrée par la figure 2. Elle est tracée pour différentes valeurs du nombre N de capteurs présents sur la zone surveillée. On constate que, plus ils sont nombreux dans une certaine limite, plus l'erreur converge rapidement. Il est à noter enfin la croissance de cette erreur à l'instant t = 100, provoquée par l'extinction de la première source et l'allumage de la seconde.

VII. Conclusion

Nous avons décrit une nouvelle classe d'algorithmes d'apprentissage non-linéaires en ligne avec contrainte de positivité sur la solution. Ceux-ci ont été appliqués au problème d'identification distribuée d'un champ scalaire positif par un réseau de capteurs. Les perspectives de ce travail concernent l'étude de la convergence de ces méthodes et, compte tenu de leur caractère générique vis-à-vis du critère optimisé, leur mise en oeuvre avec d'autres divergences.

Références

2 n

 2 dont il faut déterminer les coefficients α o = [α o 1 . . . α o L ] optimaux au sens de J. Ainsi note-t-on que le problème devient linéaire par rapport à ces derniers, permettant de reprendre le formalisme développé jusqu'ici. En particulier, pour les données poissoniennes et gaussiennes vues précédemment, les fonctions coût associées s'écrivent respectivement n, -y n, ) + y n, log y n, α κ n, n, -y n, ) 2 σ avec κ n, = [κ(x n, , x 1,1 ) κ(x n, , x 1,2 ) . . . κ(x n, , x N,L N )] conformément au Théorème de Représentation Généralisé. La contrainte de positivité de ψ o (x), pour tout x de X , se traduit alors par l'usage de noyaux positifs et la recherche d'une solution α o positive.

N 2 N

 2 n=1 (d n, -ψ ( ) (x n ))

Fig. 1 .

 1 Fig. 1. Champs de température estimés à différents instants. De gauche à droite, et de haut en bas : t = 1 (init.), t = 50, t = 100 (changement de source), t = 110, t = 150, t = 200. Les capteurs appartenant au dictionnaire sont marqués d'un point rouge, les autres d'un point bleu.

Fig. 2 .

 2 Fig.2. Evolution de l'erreur quadratique moyenne normalisée au cours du temps, pour différentes valeurs du nombre N de capteurs disponibles sur la zone surveillée.

Afin de résoudre l'équation[START_REF] Honeine | Distributed regression in sensor networks with a reduced-order kernel model[END_REF], de la forme ϕ(α) = 0, on considère le problème ϕ(α) + α = α auquel on applique une méthode de point fixe.


DISTRIBUTED LEARNING WITH KERNELS IN WIRELESS SENSOR NETWORKS FOR PHYSICAL PHENOMENA MODELING AND TRACKING

Cédric Richard (1) , Paul Honeine (2) , Hichem Snoussi (2) , André Ferrari (1) , Céline Theys (1) (1) Laboratoire Fizeau (UMR CNRS 6525, OCA), Université de Nice Sophia-Antipolis, 06108 Nice, France (1) Institut Charles Delaunay (FRE CNRS 2848), Université de technologie de Troyes,10010 Troyes, France Wireless ad-hoc sensor networks have emerged as an interesting and important research area in the last few years. They rely on sensor devices deployed in an environment to support sensing and monitoring, including temperature, humidity, motion, acoustic, etc. Low cost and miniaturization of sensors involve limited computational resources, power and communication capacities. Consequently, wireless ad-hoc sensor networks require collaborative execution of a distributed task on a large set of sensors, with reduced communication and computation burden.

In this paper, we consider the problem of modeling physical phenomena such as temperature field distributions, and track their evolution. Many approaches have been proposed in the signal processing literature to address this issue with collaborative sensor networks. In the framework of distributed inference with parametric models, focus was clearly put on determining how the channel capacity limits the quality of estimates. See, e.g., [START_REF] Swami | Wireless Sensor Networks: Signal Processing and Communications[END_REF] and references inside. Model-based techniques were also exploited to compress the data based on their temporal and spatial redundancy [START_REF] Guestrin | Distributed regression: an efficient framework for modeling sensor network data[END_REF]. Since they do not depend on some arbitrary modeling assumptions, applications of wireless sensor networks provide a strong motivation for the use of nonparametric methods for decentralized inference. See [START_REF] Predd | Distributed learning in wireless sensor networks[END_REF] for a survey. Model-independent methods based on kernel machines have recently been investigated. In particular, a distributed learning strategy has been successfully applied to regression in wireless sensor networks [START_REF]Distributed kernel regression: An algorithm for training collaboratively[END_REF]. In this paper, we take advantage of this framework to derive a new approach.

Kernel-based methods have gained wide popularity over the last decade. Initially derived for regression and classification with support vector machines, they include classical techniques such as least-squares methods and extend them to nonlinear functional approximation. We wish to determine a function ψ o (•) defined on the area of interest X that best models a distribution such as a temperature field. The latter is learned from the information coupling sensor locations and measurements. The information from the N sensors located at x n ∈ X which provide measurements d n ∈ lR, with i = n, . . . , N , is combined in the set of pairs {(x 1 , d 1 ), . . . , (x N , d N )}. Consider a reproducing kernel κ : X × X → lR. Let us denote by H its reproducing kernel Hilbert space (RKHS) with inner product • , • H . This means that every function ψ of H can be evaluated at any x ∈ X by ψ(x) = ψ, κ x H , where κ x denotes κ(•, x). This allows us to write κ(x i , x j ) = κ x i , κ x j H , which defines the reproducing property. One of the most widely used kernel is the Gaussian kernel κ(x i , x j ) = e -xi-xj 2 /2σ 2 with σ the kernel bandwidth. We seek the optimal ψ o that minimizes the following global cost function

J(ψ) = N n=1 E|ψ(x n ) -d n | 2 , that is, J(ψ) = N n=1 E| ψ, κ x n H -d n | 2 using
the reproducing property. The cost function J(ψ) can be rewritten as follows, for any node k ∈ {1, . . . , N },

J(ψ) = J k (ψ) + N n=1 n =k J n (ψ) with J k (ψ) = n∈N k c n,k E| ψ, κ x n H -d n | 2 the
c n,k = 0 if n ∈ N k , C1 1 = 1 1 and 1 1 ⊤ C = 1 1 ⊤ . Let us write ψ o n = arg min ψ∈H J n (ψ).
Finally, it is equivalent to minimize J(ψ) or the following cost function

J ℓ (ψ) = n∈N k c n,k E| ψ, κ x n H -d n | 2 + N n=1 n =k ψ -ψ o n 2 Hn
Minimizing J ℓ k (ψ) requires the nodes to have access to global information ψ o n . To facilitate distributed implementations, consider the relaxed local cost function

J r k (ψ) = n∈N k c n,k E| ψ, κ x n H -d n | 2 + n∈N k /{k} b n,k ψ -ψ o n 2 H
where ψ o n is now the best estimate available at node n. The iterative steepest-descent approach for minimizing J r k (ψ) can be expressed in the following form

ψ k,i = ψ k,i-1 - µ 2 ∇J r k (ψ k,i-1 ) with ∇J r k (ψ k,i-1 ) = n∈N k 2 c n,k E (ψ k,i-1 (x n ) -d n ) κ x n + n∈N k /{k} 2 b n,k (ψ k,i-1 -ψ o n )
Incremental algorithms are useful for minimizing sums of convex functions. They consist of iterating sequentially over each sub-gradient, in some predefined order. This leads to the following possible distributed learning strategies.

Adapt-then-Combine kernel LMS

For each time instant i and each node k, repeat

φ k,i = ψ k,i-1 -µ k n∈N k c n,k (ψ k,i-1 (x n ) -d n,i ) κ xn ψ k,i = n∈N k b n,k φ k,i

Combine-then-Adapt kernel LMS

For each time instant i and each node k, repeat

φ k,i-1 = n∈N k b n,k ψ k,i-1 ψ k,i = φ k,i-1 -µ k n∈N k c n,k (φ k,i-1 (x n ) -d n,i ) κ xn
∂Θ(x, t) ∂t -c∇ 2 x Θ(x, t) = Q(x, t).
Here Θ(x, t) denotes the temperature as a function of space and time, c is a medium-specific parameter, ∇ 2 x is the Laplace spatial operator, and Q(x, t) is the heat added. We studied the question of monitoring the evolution of the temperature field distribution in a square region with open boundaries and conductivity c = 0.1, using N = 100 sensors deployed randomly on a grid. Given some measurements d n,i = Θ(x n , t i ) + z n,i , with z n,i a i.i.d. Gaussian noise, the problem was to estimate the temperature Θ(x n , t i ) via ψ i (x n ) based on the Gaussian kernel.

Two heat sources of intensity 200 W were placed within the region, the first one was activated from t = 1 to t = 100, and the second one from t = 100 to t = 200. The experimental setup will be described in the camera-ready paper. Fig. 1 illustrates the estimated temperature field distribution at different times by the Adapt-then-Combine strategy, for several values of |N k |. The convergence of the proposed algorithm is illustrated in Fig. 2 where we show the evolution over time of the normalized mean-square prediction error. The abrupt change in heat sources at t = 100 is clearly visible, and highlights the convergence behavior of the proposed algorithm. Comparison with an iterative distributed gradient algorithm clearly shows a faster convergence speed.

  local cost function, and N k the set of neighbors of sensor k. Here, c n,k is the (n, k)-th component of C, which satisfies the following constraints:

Fig. 1 .

 1 Fig. 1. Snapshots of the evolution of the estimated temperature at t = 100 (left), t = 150 (center) and t = 200 (right). Sensors are shownwith small blue dots. Big red dots represent the sensors of interest according to a criterion that will be described in the camera-ready paper.

4 Fig. 2 .

 42 Fig. 2. Learning curve obtained from t = 1 to t = 200. Time t = 100 corresponds to a system modification.


Jie Chen

Cédric Richard

Paul Honeine

José Carlos

M Bermudez

Non-negative Distributed Regression for Data Inference in Wireless Sensor Networks

Wireless sensor networks are designed to perform on inferences about the environment they are sensing. Due to the inherent physical characteristics of systems under investigation, non-negativity is a desired constraint that can be imposed on the system parameters in some real-life phenomena sensing tasks. In this paper, we propose a kernel-based machine learning strategy to deal with regression problems. Multiplicative update rules are derived in this context to ensure the non-negativity constraints to be satisfied. A distributed algorithm that requires only communication between neighbors is proposed to cope with typical limited energy and bandwidth resources. Synthetic data managed by heat diffusion equations are used to compare the proposed and known algorithms and to illustrate their tracking capabilities.

I. INTRODUCTION

Wireless sensor networks (WSNs) rely on sensor devices deployed in an environment to provide an inexpensive way to monitor physical phenomena. In many real-life phenomena, including biological and physical ones, physical characteristics inherent to the system under investigation require the imposition of non-negativity constraints on the parameters to estimate. For instance, observations in studies of concentration fields or thermal radiation fields are always described with non-negative values ( in ppm or in Kelvin). Non-negativity as a physical constraint has received growing attention from the signal processing community during the last decade [START_REF] Saul | Statistical signal processing with nonnegativity constraints[END_REF].

Non-parametric approaches based on reproducing kernel methods have recently been successfully applied to distributed regression with collaborative networks. In [START_REF] Guestrin | Distributed regression: an efficient framework for modeling sensor network data[END_REF], the authors present a general framework for distributed linear regression motivated by WSNs. In [START_REF] Predd | Regression in sensor networks: Training distributively with alternating projections[END_REF], a learning algorithm based on successive orthogonal projections is derived to solve the regularized kernel least-squares problem for regression in sensor networks. In [START_REF] Honeine | Functional estimation in Hilbert space for distributed learning in wireless sensor networks[END_REF], the authors present a projection based kernel distributed learning strategy with reduced order models by using a sparsification criterion. These works provide comprehensive studies in the functional regression and estimation for distributed learning in WSNs. In [START_REF] Chen | System identification under non-negativity constraints[END_REF], the authors explore online algorithm for linear system identification. However, none of these algorithms can be used directly to solve the estimation problems in sensor networks under non-negativity constraints.

In this paper, we concentrate on the problem of modeling physical phenomena under non-negativity constraints, and of tracking its evolution. Firstly we formulate the non-negative regression with kernels in a centralized context. A simple multiplicative algorithm is derived to solve this problem. Then we show how the optimization problem can be relaxed to a problem of distributed regression in which nodes only need to communicate with neighbors.

II. NON-NEGATIVE REGRESSION FOR INFERENCE

Within the context of learning in a wireless sensor network of N sensors, we often model a physical phenomenon as a function of the location. Consider a relationship ψ(.) between the sensor's measurement and its position x n . We seek to estimate the function ψ(•) based on newly available positionmeasurement data y n to minimize the summed square error

min ψ∈H N ∑ n=1 E (ψ(x n ) -y n ) 2 . ( 1
)
By virtue of the representer theorem, the function ψ(•) of reproducing Hilbert kernel space H can be written with a kernel expansion ψ(•) = ∑ N n=j α j κ(•, x n ). Doing that, the cost function can be written as

J(α) = N ∑ n=1 E( N ∑ j=1 α j κ(x n , x j ) -y n ) 2 = N ∑ n=1 E ( α ⊤ κ xn -y n ) 2 . ( 2
)
After determining the weight vector α , the field can be inferred at any points x. One of the most widely used kernels is the Gaussian kernel κ(x i , x j ) = e -∥xi-xj ∥ 2 /2σ 2 When a non-negative field is to be estimated, and considering that the Gaussian kernel is always positive, each component of the coefficient vector α should be constrained to be nonnegative to ensure a non-negative inference function ψ(x) at any given position x. The constrained optimization problem can be formalized as

α o = arg min α J(α) (3) subject to α ≥ 0 (4)
The gradient of J(α) is easily computed as follows

∇J(α) = N ∑ n=1 E ( κ xn κ ⊤ xn α -y n κ xn ) (5)
As the evaluation of the gradient usually cannot be achieved in many real-life applications, we use the instantaneous estimator

∇J(α) = N ∑ n=1 ( κ ⊤ xn κ xn α -y n κ xn ) (6)
Note that ψ(.) is linear with respect to the kernel functions κ(., x n ), although it is nonlinear with respect to x n . Therefore, we use an algorithm similar to that proposed in our previous work [START_REF] Chen | System identification under non-negativity constraints[END_REF] to solve the problem. We decompose the gradient -∇J(α) as follows,

[-∇ J(α(k))] i = [U (α(k))] i -[V (α(k))] i (7)
where [U (α(k))] i and [V (α(k))] i are strictly positive components. The update equation for the i-th component can then be expressed as

α i (k + 1) = α i (k) [U (α(k))] i [V (α(k))] i ( 8
)
This expression is referred to as the multiplicative weight update algorithm. If we initialize the weight vector with a positive vector, the constraints will be always satisfied due to the non-negativity of [U (α(k))] i and [V (α(k))] i . The gradient defined by (6) can be decomposed as in (7) by setting

U (α(k)) = ∑ N n=1 y n κ xn + ξ and V (α(k)) = ∑ N n=1 κ κ ⊤ xn α(k) + ξ, with ξ positive to avoid [U (α(k))
] i to become negative due to the noise. The centralized algorithm vector weight update is then

α(k + 1) = α(k)diag ( [ ∑ N n=1 y n κ xn + ξ] i [ ∑ N n=1 κ xn κ ⊤ xn α(k) + ξ] i ) (9)

III. DISTRIBUTED REGRESSION WITH DIFFUSION STRATEGY IN WSNS

In what follows, we show how the optimization problem in (2) can be relaxed for distributed inference. Let N denote the set of neighbors for sensor k. Consider an N × N matrix B with entries {b n,k } defined by b n,k = 0 if n / ∈ N k and B1 = 1, B = 1 ⊤ . With the constraint of communication range, the cost function of ( 2) is rewritten as follows

J(α) = J k (α) + N ∑ n=1,n̸ =k J n (α) (10)
We define diagonal matrices C k for each node k with elements c k,i,i = 1 if i ∈ N k , and c k,i,i = 0 otherwise. The local cost function J k (α) is defined as

J k (α) = N ∑ n=1 b n,k E ( α ⊤ C n κ xn -y n ) 2 (11)
which is actually equivalent to

J k (α k ) = ∑ n∈N k b n,k E ( α ⊤ k κ xn -y n ) 2 (12)
Each node could only communicate with the nodes within the range of a neighborhood of node k. The instananeous gradient of the cost function (10) at each node k is

[∇J(α k )] i = [ ∑ n∈N k b n,k ( κ xn κ ⊤ xn α k -y n κ xn ) ] i +   N ∑ n=1,n̸ =k C n ∇J n (α)   i ( 13
)
where i ∈ N k . To relax the problem so that sensors only need to get information from its neighbors, we use

[∇J(α k )] i = [ ∑ n∈N k b n,k ( κ xn κ ⊤ xn α k -y n κ xn ) ] i +   N ∑ n∈N k ,n̸ =k C n ∇J n (α)   i ( 14
)
The first part of ( 14) can be viewed as the gradient of local cost function ∇J k (α k ). Using the proposed multiplicative algorithm,

[-∇J k (α k )] i is decomposed into two positive components [U k (α(k))] i = [∑ n∈N k b n,k y n κ xn ] i + ξ and [V k (α(k))] i = [∑ n∈N k b n,k κ xn κ ⊤ xn α k (k) ] i + ξ.
The second part of ( 14) can be viewed as a regularization item which is decomposed using

[ U k (α(k)) ] i = ∑ n∈N k ,n̸ =k [U n (α(k))] i and [ V k (α(k)) ] i = ∑ n∈N k ,n̸ =k [V n (α(k))] i , where [U (α(k))] i and [V (α(k))] i are transferred from its neighbors.
Finally, the coefficients update rule for node k is written in multiplicative form as

α i (k + 1) = α i (k) [U k (α(k))] i + [ U k (α(k)) ] i [V k (α(k))] i + [ V k (α(k)) ] i (15)
IV. SIMULATION EXPERIMENTS We consider a classical application of estimating a heat diffusion field governed by the partial differential equation

∂T (x, t) ∂t -c∇ 2 x T (x, t) = Q(x, t).
We studied the problem of monitoring the evolution of the heat using N = 100 random sensors. One heat source was activated from t = 1 to t = 100, and the second one from t = 100 to t = 200. The bandwidth of Gaussian function σ = 0.1826. Two random scenarios, respectively depicted in Figure 1 and Figure 2, are taken into simulations. In the 1 st scenario, two sources are located at positions relatively "poor" of sensors; otherwise, in the 2 nd scenario, the two sources are located at "better" positions. The algorithms compared are: 1) Centralized multiplicative algorithm; 2) Proposed distributed multiplicative; 3) Centralized Barzilai-Borwein gradient projection; 4) Distributed gradient projection. The learning curves are shown in figures 3 and 4 respectively. The abrupt change in heat sources at t = 100 is clearly visible, and the proposed algorithm has good performance with little computational burden.

V. CONCLUSION

In many real-life phenomena non-negativity is a desired constraint that must be imposed on the parameters to estimate due to the inherent physical characteristics of systems In this paper, we proposed a multiplicative method for data inference under non-negativity constraints. Under the context of wireless sensor networks, we developed a distributed learning algorithm to enable each sensor to estimate the non-negative field with the help of neighbor information. The proposed algorithm also shows a good performance in its tracking capacity.

Fig. 1 .

 1 Fig. 1. 1st simulation scenario.The edges between two nodes show the neighborhood relation. Two magenta △ represent the positions of two sources.

Fig. 2 .

 2 Fig.2. 2nd simulation scenario. In this scenario, the positions of two sources are more close to sensors than that in scenario1.

Fig. 3 .

 3 Fig. 3. Convergence comparison of the 1st scenario.

Fig. 4 .

 4 Fig. 4. Convergence comparison of the 2nd scenario.


Paul Honeine
email: paul.honeine@utt.fr

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

A SIMPLE SCHEME FOR UNMIXING HYPERSPECTRAL DATA BASED ON THE GEOMETRY OF THE N-DIMENSIONAL SIMPLEX

In this paper, we study the problem of decomposing spectra in hyperspectral data into the sum of pure spectra, or endmembers. We propose to jointly extract the endmembers and estimate the corresponding fractions, or abundances. For this purpose, we show that these abundances can be easily computed using volume of simplices, from the same information used in the classical N-Findr algorithm. This results into a simple scheme for unmixing hyperspectral data, with low computational complexity. Experimental results show the efficiency of the proposed method.

INTRODUCTION

Based on observation satellites, remote-sensing with hyperspectral image analysis has gained wide popularity, with applications to mineral exploration, environmental monitoring, and military surveillance. With a large number of measured wavelength bands, each pixel has a complete spectrum. It is usually assumed that each spectrum can be linearly decomposed into spectra, the so-called endmembers, each provided by some pure material. For a given image, the hyperspectral unmixing problem involves, on the one hand the extraction of these pure spectra and, on the other the estimation of their abundances within each pixel. The number of endmembers is often assumed known; otherwise methods such as the one presented in [START_REF] Chang | Estimation of number of spectrally distinct signal sources in hyperspectral imagery[END_REF] can be used to estimate it.

The first step is the endmember extraction problem. From convex geometry, this problem is equivalent to finding a dataenclosing simplex, as proposed in [START_REF] Boardman | Automating spectral unmixing of aviris data using convex geometry concepts[END_REF][START_REF] Winter | N-findr: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF]. One of the most popular (automated) algorithms for endmember extraction is Winter's N-Findr [START_REF] Winter | N-findr: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF]. The main driving force behind this algorithm is a simple iterative and scalable scheme. With a random initial set defining the simplex, it seeks to stretch the simplex in order to maximize its volume by visiting each pixel. The N-Findr is a fast endmember extraction technique with low computational cost, making it the most widely used method despite the restriction on belonging endmembers in the image. The second step involves estimating the abundance of each endmember for all pixels. Many techniques have been proposed to estimate the fraction of pure spectra in a given spectrum. The main difficulties reside in constraints on the estimated coefficients in order to represent the physical abundancies, with the sum-to-one and the non-negativity constraints(see [START_REF] Keshava | An algorithm taxonomy for hyperspectral unmixing[END_REF][START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF] and references therein).

In this paper, we give a simple algorithm for extracting the endmembers and their respective adundances, jointly in a single step. The proposed approach takes advantage of the well-known N-Findr scheme for endmember extraction. In the same spirit, we show that we can provide an efficient iterative algorithm for estimating the abundancies. Both the endmember extraction and the abundancies estimation are done in a single step, using the concept of barycentric coordinates which are computed from volumes of simplices. This is done without any additional computational cost as opposed to classical least square estimation techniques [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF].

The linear (un-)mixing model

In a hyperspectral image, the hypothesis of a linear mixed model is often used. The spectrum of a pixel is given as a linear combination of some pure spectra, called endmembers. Let x ℓ be the spectrum of the ℓ-th pixel, then the mathematical model is

x ℓ = n i=1 α i x ωi + ǫ, (1)
where endmembers are defined by spectra x ωi , assumed represented by some pixels {ω 1 , ω 2 , . . . , ω n }, and ǫ corresponds to the unfitness of the model, often treated as a Gaussian noise. In order for the coefficients α 1 , α 2 , . . . , α n to represent the physical abundance fraction associated to each endmember, two constraints must be imposed on theses coefficients:

• The sum-to-one constraint, called hereafter the equality constraint, with

n i=1 α i = 1.
• The non-negativity constraint, with

α 1 , α 2 , . . . , α n ≥ 0.
Estimating the coefficients based on the model (1) subject to both constraints requires advanced optimization techniques, as studied for instance in [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF]. In this paper, we give a direct scheme to solve the equality-constrained optimization problem, incorporating naturally the endmember extraction procedure. Furthermore, we provide a geometric interpretation to the violation of the non-negativity constraint.

THE METHOD

Unmixing hyperspectral data based on the linear model ( 1) involves two tasks: extracting the endmembers x ωi and computing their coefficients α i , for i = 1, 2, . . . , n, and for the spectrum of each pixel. We begin next with the extraction method, as given in the N-findr scheme [START_REF] Winter | N-findr: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF], and then derive in next paragraph the method to estimation the corresponding coefficients.

But before, we give a definition of the volume of a simplex. Let X = {x ω1 , x ω2 , . . . , x ωn } be the set of estimated endmembers, with ω i ∈ {1, 2, . . .}. The oriented volume of the simplex defined by the vertices X is given by

V X = 1 (n -1)! det 1 1 • • • 1 x ω1 x ω2 • • • x ωn , ( 2
)
where det is the determinant operator. While taking the absolute value of this expression gives the (conventional, positivevalued) volume of a simplex, the virtues of the oriented (signed) volume will be demonstrated in estimating the coefficients.

Endmember extraction

We take advantage of the fact that the endmembers define the vertices of a simplex englobing all spectra of the image. Thus, this simplex has the largest volume among all simplices constructed from other spectra. Therefore, one seeks the simplex of largest volume, in an iterative manner by visiting each pixel. At initialization, random pixels are selected as endmembers 1 . The following process is iterated for each pixel, where x ℓ is its spectrum: One at a time, each endmember is replaced by the pixel under investigation, and the oriented volume of the resulting simplex is computed. Let

V X\{xω i }∪{x ℓ }
denotes the simplex (oriented) volume of vertices X\{x ωi }∪ {x ℓ }, with \ the set difference defined by the set X with x ωi removed and x ℓ added. Besides the initial set with simplex volume V X , we have n candidate sets, each defining a simplex volume V X\{xω i }∪{x ℓ } . By comparing these volumes, two cases can be distinguished:

1 One may also operate an appropriate selection of the starting set, as proposed for instance in [START_REF] Plaza | An improved n-findr algorithm in implementation[END_REF]. However, such improvements are beyond the scope of this communication.

• if max i |V X\{xω i }∪{x ℓ } | < |V X |, then the initial set of endmembers remains unchanged;

• otherwise, an entry of the initial set is substituted with x ℓ to give the new endmembers set. The outgoing spectrum ω i is identified such as

i = arg max |V X\{xω i }∪{x ℓ } |.

Coefficients estimation

We are now in a position to compute the coefficients. This is done at the same iteration, using the computed oriented volumes. Given the simplex defined by the vertices x ω1 , x ω2 , . . . , x ωn , any x ℓ can be written as a linear combination of these vertices, i.e.

x ℓ = n i=1 α i x ωi , (3)
subject to

n i=1 α i = 1 (4)
To solving this constrained optimization problem, we combine the above expressions into the following augmented2 linear system

1 1 • • • 1 x ω1 x ω2 • • • x ωn α = 1 x ℓ , where α = [α 1 α 2 • • • α n ] ⊤
is the vector of coefficients. Using the well-known Cramer's rule, the solution of this linear system can be expressed in terms of the determinants of the above matrix and of matrices obtained from it with one column substituted by the right-hand-side vector. Therefore, we can write

α 1 = det 1 1 • • • 1 x ℓ x ω2 • • • x ωn det 1 1 • • • 1 x ω1 x ω2 • • • x ωn , α 2 = det 1 1 • • • 1 x ω1 x ℓ • • • x ωn det 1 1 • • • 1 x ω1 x ω2 • • • x ωn , . . . α n = det 1 1 • • • 1 x ω1 x ω2 • • • x ℓ det 1 1 • • • 1 x ω1 x ω2 • • • x ωn . x ℓ x ω2 x ω1 x ω3 x ℓ x ω2 x ω1 x ω3
Fig. 1. Illustration of the simplex (here triangle) in a 2-dimensional Euclidean space. Left: the current data x ℓ is inside the simplex defined by {xω 1 , xω 1 , xω 1 }, and thus all coefficients are non-negative. Right: the data is outside the simplex and thus has at least one negative coefficient. Thus, each coefficient α i equals a ratio of the oriented volumes of two simplices, the one with vertices X given by endmembers and the one resulting from the latter by substituting vertex x ωi with x ℓ . In other words, we have

α i = V X\{xω i }∪{x ℓ } V X , (5)
for all i = 1, 2, . . . , n. Since these volumes are initially computed in order to stretch the initial simplex for enclosing most spectra, no additional computational cost is required to evaluate the coefficients at each iteration.

On the non-negativity constraint

From the literature of geometry of convex sets, the coefficients computed using (5) are called (homogeneous) barycentric coordinates or areal coordinates, defined by the vertices X of the simplex. It is well known that x ℓ is inside the simplex if and only if all α i 's are non-negative; otherwise, when it is outside, there exists at least one negative coefficient. The corresponding simplices are illustrated in Fig. 1 for two-dimensional data in two case, whether x ℓ lies inside or outside the endmembers simplex. Therefore, when the solution of the constrained optimization problem (3)-( 4) gives at least one negative α i , this means that the x ℓ is outside the simplex, and thus cannot be written in terms of linear combination of its vertices with both equality and non-negativity constraints. It is worth noting that such limitation is valid for any simplex-based approach. Therefore, this provides a measure of adequacy of the geometric approach, i.e. the simplex englobing the data with the simplex vertices available within the data.

EXPERIMENTATIONS

In order to illustrate the proposed method, we simulated a synthetic hyperspectral image from a linear combination of three pure materials. These materials are available from the USGS library [START_REF] Clark | USGS digital spectral library splib06[END_REF], and correspond to golden grass, red brick, and cedar, with spectra illustrated in Fig. 2 (left). A 128-by-128 hyperspectral image is generated from pixels given by the model [START_REF] Chang | Estimation of number of spectrally distinct signal sources in hyperspectral imagery[END_REF], where ǫ corresponds to a white Gaussian noise of variance 0.01. Using the same noisy model with the endmembers, they are incorporated in the image with the canonical coefficients {(1, 0, 0); (0, 1, 0); (0, 0, 1)}, as illustrated in Fig. 2

(right).

Most simplex-based methods require a preprocessing dimensionality reduction technique. For this purpose, we apply a classical PCA as suggested in [START_REF] Winter | N-findr: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF]; alternatively, more dedicated methods such as minimum noise fraction technique can be considered. We apply the proposed method with the number of endmembers set to 3. The resulting largest-volume simplex encloses 71% of the spectra, following from incorporating noise, not only in the coefficients by also in the endmembers. To measure the performance of the algorithm, we compute the spectral angle error [START_REF] Kruse | The spectral image processing system (SIPS): Interactive visualization and analysis of imaging spectrometer data[END_REF] for each pixel, between its initial spectrum and the one computed from the obtained coefficients, of the form α(x i , x j ) = cos -1

x i , x j x i x j .

The histogram of these errors given in Fig. 3 shows a small angular error for the hyperspectral image.

CONCLUSION

In this paper, we studied the problem of unmixing hyperspectral data. Based on the well-known Winter's N-Findr scheme, we showed that we can jointly operate extraction of endmembers and computing the abundances, using volumes of simplices.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Endmember spectra: (left) initial spectra and (right) their noisy versions used in generating the hyperspectral image, including two random mixed spectra.

Using an augmented system has been previously introduced in the literature, with the fixed-point-free transform[START_REF] Craig | Minimum-volume transforms for remotely sensed data[END_REF]. To our knowledge, this is the first time that barycentric coordinates are applied to provide a simple optimization scheme.


Paul Honeine

José Carlos M Bermudez

Hichem Snoussi

Mehdi Essoloh

FUNCTIONAL ESTIMATION IN HILBERT SPACE FOR DISTRIBUTED LEARNING IN WIRELESS SENSOR NETWORKS

Keywords: Intelligent sensors, adaptive estimation, distributed algorithms, nonlinear systems

In this paper, we propose a distributed learning strategy in wireless sensor networks. Taking advantage of recent developments on kernel-based machine learning, we consider a new sparsification criterion for online learning. As opposed to previously derived criteria, it is based on the estimated error and is therefore is well suited for tracking the evolution of systems over time. We also derive a gradient descent algorithm, and we demonstrate its relevance to estimate the dynamic evolution of temperature in a given region.

INTRODUCTION

Wireless ad-hoc sensor networks have emerged as an interesting and important research area in the last few years. They rely on sensor devices deployed in an environment to support sensing and monitoring, including temperature, humidity, motion, acoustic, etc. Low cost and miniaturization of sensors involve limited computational resources, power and communication capacities. Consequently, wireless adhoc sensor networks require collaborative execution of a distributed task on a large set of sensors, with reduced communication and computation burden.

In this paper, we consider the problem of modeling physical phenomena, such as a temperature field, and track its evolution. Many approaches have been proposed in the signal processing literature to address this issue with collaborative sensor networks. See [START_REF] Predd | Distributed learning in wireless sensor networks[END_REF] for a survey. As explained in [START_REF] Predd | Distributed learning in wireless sensor networks[END_REF], the incremental subgradient optimization scheme derived in [START_REF] Rabbat | Distributed optimization in sensor networks[END_REF] for (a single) parameter estimation is not appropriate for large-order models. In [START_REF] Guestrin | Distributed regression: an efficient framework for modeling sensor network data[END_REF], the authors use both spatial correlation and time evolution of sensors to propose a reduced-order model. However, this approach highly depends on the modeling assumption. Recently, model-independent methods have been investigated. A distributed learning strategy in sensor networks is studied in [START_REF] Predd | Distributed kernel regression: An algorithm for training collaboratively[END_REF], where each sensor acquires information from neighboring sensors to solve locally a least-squares problem. Unfortunately, this broadcast leads to high energy consumption.

Recently, kernel machines for nonlinear functional learning have gained popularity [START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF][START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF]. Nevertheless, these methods are not suitable for distributed learning in sensor networks as the order of models scales linearly with the number of deployed sensors and measurements. In order to circumvent this drawback, we propose in this paper to design reduced order models by using an easy to compute sparsification criterion. As opposed to a criterion previously derived 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][START_REF] Honeine | Distributed regression in sensor networks with a reduced-order kernel model[END_REF], it depends on the estimated error. This approach is, therefore, more relevant in updating the model since it is based on available measurements. Based on this criterion and a projection scheme, we derive the learning algorithm by incrementing the model order if necessary, leaving it unchanged, or even decreasing it. We illustrate the proposed approach for learning a temperature field and tracking its evolution over time. Before proceeding, we briefly review functional learning with kernels and its online setting.

ONLINE LEARNING WITH KERNELS

Consider a reproducing kernel κ : X × X → IR. Let us denote by H its reproducing kernel Hilbert space (RKHS) with inner product • , • H . This means that every ψ(•) of H can be evaluated at any x ∈ X by ψ(x) = ψ(•), κ(•, x) H . This allows us to write κ(x i , x j ) = κ(•, x i ), κ(•, x j ) H , which defines the so-called reproducing property. One of the most widely used reproducing kernel is the Gaussian kernel, given by κ(x i , x j ) = e -xi-xj 2 /2σ 2 with σ the kernel bandwidth.

Within the context of distributed learning in a wireless sensor network, we model a physical phenomenon, e.g., a temperature field, as a function of the location x. Let us denote it by ψ n (•) ∈ H where X represents the 2-D space. We seek to estimate the function ψ n (•) at sensor n based on newly available position-measurement data, (x n , d n ), and the previous estimate ψ n-1 (•). For this purpose, we consider the following problem

ψ n = arg min ψ∈H ψ n-1 -ψ 2 H (1) subject to ψ n (x n ) = d n .
(

This optimization problem can be interpreted as a classical adaptive filtering problem, applied here to functional estimation in a RKHS. Expression (1) corresponds to the classical principle of minimum disturbance, and the constraint (2) sets to zero the a posteriori error. Though a large class of adaptive filtering techniques can be used here, we restrict ourselves to a gradient descent approach as studied in [START_REF] Smale | Online learning algorithms[END_REF] and we consider the updating step

ψ n = ψ n-1 + η n (d n -ψ n-1 (x n ))κ(x n , •).
In what follows, we set the tunable positive stepsize to η n = 1 as used in [START_REF] Dodd | Function estimation in Hilbert space using sequential projections[END_REF]. In addition, we consider unit-norm kernel functions, i.e., κ(x, x) = 1 for any x ∈ X . The above expression yields the updating rule

ψ n = ψ n-1 + ǫ n κ(x n , •), (3)
where

ǫ n = d n -ψ n-1 (x n ) is the a priori estimation error.
Applying this updating rule sequentially to n sensors leads to the n-order model

ψ n = n i=1 α i κ(x i , •), (4)
where all the coefficients α i are identical to those of ψ n-1 , except α n = ǫ n .

The above updating rule leads to models with orders equal to the number of available data. Such models are not suitable for large-scale data problems or online learning. This is the case of models generated by most kernel machines. To overcome this drawback, we propose in the next section a new online sparsification technique to control the model order.

THE PROPOSED SPARSIFICATION CRITERION

We consider an m-order model, with m several orders of magnitude lower than n, defined by

ψ n (•) = m k=1 α k κ(x ω k , •), (5)
where {ω 1 , . . . , ω m } is a subset of {1, . . . , n}. We thus restrict the expansion to m kernel functions carefully selected among the n kernel functions in model [START_REF] Predd | Distributed kernel regression: An algorithm for training collaboratively[END_REF].

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 proposed a sparsification technique for designing models with kernel functions having small coherence, the latter being defined as

max i =j | κ(x ωi , •), κ(x ωj , •) H |/ κ(x ωi , •) H κ(x ωj , •) H .
The sparsification rule consisted of including, for each sensor n, the kernel function κ

(x n , •) into the model if max k=1,...,m |κ(x n , x ω k )| κ(x n , x n )κ(x ω k , x ω k ) ≤ ν 0 , (6)
with ν 0 a threshold in [0, 1[ determining the level of sparsity of the model. 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 studied this sparsification rule for online learning. We also derived some properties of the resulting model as well as connections to other sparsification techniques. In [START_REF] Honeine | Distributed regression in sensor networks with a reduced-order kernel model[END_REF], we investigated the application of this criterion for wireless sensor networks. The independence of the sparsification rule with respect to measurements and to estimation errors limited the performance of the resulting function estimation process for that application. In this paper, we propose to overcome this limitation by using the concept of coherence between the ψ k 's.

The function ψ n defined in (3) is selected as the new model if max k=1,...,n-1

| ψ n , ψ k H | ψ n H ψ k H ≤ ν, (7)
with ν a threshold. Otherwise, we use the projection of ψ n onto the space H m-1 spanned by the m -1 previously added kernel functions. It is obvious that solving this problem is untractable in practice since we need to know all previous estimated functions, ψ 1 , ψ 2 , . . . , ψ n-1 . However, because these functions belong to H m-1 , we can circumvent this difficulty as explained below.

Proposition 1. Let ψ ⊥ n be the projection of ψ n onto the space spanned by the m -1 kernel functions. If we have

ψ n , ψ ⊥ n H ψ n H ψ ⊥ n H ≤ ν, (8)
then the inequality (7) is satisfied.

Sketch of proof. To prove this, note that

ψ ⊥ n = arg max φ∈Hm-1 ψ n , φ H ψ n H φ H .
Since the estimated functions ψ 1 , ψ 2 , . . . , ψ n-1 belong to the space H m-1 , the criterion (8) directly leads to [START_REF] Honeine | On-line nonlinear sparse approximation of functions[END_REF].

Upon the arrival of a new data (x n , d n ), one of the following two alternatives holds. If ( 8) is satisfied, the kernel function κ(x n , •) is then added to the model according to (3). Otherwise, the model order is not incremented and we consider the closest function to ψ n in H m-1 , that is, ψ ⊥ n . Additionally to this rule, we propose a strategy to decrease the model order. With sensors being revisited in order to follow the evolution of the system over time, new data may correspond to a sensor1 that was incorporated in the model in a previous pass.

Let κ(x n , •) be a kernel function that is already in the model. Its relevance depends now on the new measurement d n . In that case, criterion ( 8) is evaluated to determine whether this kernel function should be kept or removed from the model.

According to [START_REF] Guestrin | Distributed regression: an efficient framework for modeling sensor network data[END_REF], it clearly appears that this rule depends on the estimated error. It is thus related to d n as opposed to rule [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF]. It can be shown that the order of the model resulting from rule (8) remains finite as n goes to infinity, even when the decreasing scheme is not used. Due to limited space, the proof of this property is beyond the scope of this paper.

ONLINE LEARNING ALGORITHM

In this section, we derive our online learning algorithm, with recursive techniques for both incremental and decremental stages. Before proceeding, we formulate the projection problem in a RKHS.

Projection in a RKHS

Let ψ ⊥ n = m-1 i=1 β i κ(x ωi ,
•) be the projection of ψ n defined by equation ( 3) onto the space spanned by the (m -1) kernel functions κ(x ω1 , •), . . . , κ(x ωm-1 , •). The function

ψ ⊥ n is ob- tained by minimizing ψ n -ψ ⊥ n 2
H with respect to the β i 's, namely,

ǫ n κ(x n , •) - m-1 i=1 (β i -α i ) κ(x ωi , •) 2 H .
By expressing this norm in terms of inner products and using the reproducing property, we formulate the optimization problem as

min β (β -α) ⊤ K m-1 (β -α) + ǫ 2 n -2ǫ n (β -α) ⊤ κ n ,
where α, β and κ n are (m -1)-length column vectors with entries α i , β i , and κ(x ωi , x n ), respectively, and K m-1 is a (m -1)-by-(m -1) matrix whose (i, j)-th entry is given by κ(x ωi , x ωj ). By taking the derivative of the above objective function with respect to β, and setting it to zero, we get

β = α + ǫ n K -1 m-1 κ n , (9)
where we have assumed that the Gram matrix K m-1 is nonsingular. We can now present the different building blocks of the algorithm.

The sparsification criterion

The sparsification criterion needs to be evaluated by each sensor node n. The corresponding kernel function κ(x n , •) is added to the model if it satisfies the rule [START_REF] Richard | Online prediction of time series data with kernels[END_REF]. If it already belongs to the model, this rule is used to verify whether it can be removed or not. By expanding each term in the left-hand side of expression (8), we get the rule

α ⊤ K m-1 α + 2ǫ n α ⊤ κ n + ǫ 2 n κ ⊤ n K -1 m-1 κ n α ⊤ K m-1 α + 2ǫ n α ⊤ κ n + ǫ 2 n ≤ ν 2
This expression as well as equation ( 9) require to compute the inverse of the Gram matrix K m-1 . This operation can be performed by using a rank-one update, which requires O(m 2 ) operations, as derived next for both incremental and decremental stages.

Incremental and decremental steps

Increasing the model order by including κ(x n , •) into the kernel expansion requires augmenting the Gram matrix as follows

K m = K m-1 κ n κ n ⊤ κ(x n , x n ) , (10)
with κ(x n , x n ) = 1. The inverse of K m can be computed by using the rank-one update given by

A B C D -1 = A -1 0 0 0 + -A -1 B I × (D -CA -1 B) -1 -CA -1 I , (11)
with I the identity matrix. We obtain the updating rule

K -1 m = K -1 m-1 0 m-1 0 m-1 ⊤ 0 + 1 1 -κ n ⊤ K -1 m-1 κ n × -K -1 m-1 κ n 1 -κ n ⊤ K -1 m-1 1 ,
where 0 m-1 is a (m -1)-length column vector of zeros.

In the decremental stage, κ(x n , •) is removed from the model. This reduces the model order from m to m -1. The Gram matrix K m-1 is obtained from K m by considering expression [START_REF] Smale | Online learning algorithms[END_REF], where the latter matrix is arranged in order that its last column and row have entries relative to x n . Using the notation

K -1 m = Q m-1 q q ⊤ q 0 ,
we obtain from [START_REF] Dodd | Function estimation in Hilbert space using sequential projections[END_REF] the following matrix update equation

K -1 m-1 = Q m-1 - q q ⊤ q 0 .

SIMULATION RESULTS

To illustrate the relevance of the proposed technique, we consider a classical application of estimating a temperature field governed by the partial differential equation 2

∂T (x, t) ∂t -c∇ 2 x T (x, t) = Q(x, t).

2 Data simulated using MATLAB's PDE toolbox. Here T (x, t) denotes the temperature as a function of space and time, c is a medium-specific parameter, ∇ 2

x is the Laplace spatial operator, and Q(x, t) is the heat added. We studied the problem of monitoring the evolution of the temperature in a 2-by-2 square region with open boundaries and conductivity c = 0.1, using N = 100 sensors deployed randomly on a grid. Two heat sources of intensity 200 W were placed within the region, the first one was activated from t = 1 to t = 100, and the second one from t = 100 to t = 200.

Preliminary experiments were conducted to tune the parameters, yielding σ = 0.5 and ν = 0.995. In order to refine the results, 10 passes through the network were conducted at each instant t. Fig. 1 illustrates the estimated temperature field at different times. It is can be observed that the selected sensors for each snapshot follows the dynamic behavior of the heat sources. The convergence of the proposed algorithm is illustrated in Fig. 2 where we show the evolution over time of the normalized mean-square prediction error, defined on all the sensors by

1 N N n=1 (d n -ψ n-1 (x n )) 2 d 2 n .
The abrupt change in heat sources at t = 100 is clearly visible, and highlights the convergence behavior of the algorithm.

CONCLUSION

In this paper, we proposed an online learning algorithm for wireless sensor networks. It consisted of a kernel machine associated with a new sparsification criterion. We highlighted the relevance of this criterion and derived a learning algorithm with model-order control. Applications to temperature tracking with dynamic heat sources were considered, and simulation results showed the relevance of the proposed approach.

Fig. 1 .

 1 Fig. 1. Snapshots of the evolution of the estimated temperature at t = 100 (left), t = 150 (center) and t = 200 (right). Selected sensors at these instances are shown with big red dots, whereas the remaining sensors are represented by small blue dots.

3 Fig. 2 .

 32 Fig. 2. Learning curve obtained from t = 1 to t = 200. Time t = 100 corresponds to a system modification.

Sensors are assumed motionless; Otherwise, one may include a tolerance range for the positions. However, this is beyond the scope of this paper.


Paul Honeine

Cédric Richard

Mehdi Essoloh

Hichem Snoussi

LOCALIZATION IN SENSOR NETWORKS -A MATRIX REGRESSION APPROACH

In this paper, we propose a new approach to sensor localization problems, based on recent developments in machine leaning. The main idea behind it is to consider a matrix regression method between the ranging matrix and the matrix of inner products between positions of sensors, in order to complete the latter. Once we have learnt this regression from information between sensors of known positions (beacons), we apply it to sensors of unknown positions. Retrieving the estimated positions of the latter can be done by solving a linear system. We propose a distributed algorithm, where each sensor positions itself with information available from its nearby beacons. The proposed method is validated by experimentations.

INTRODUCTION

In ad-hoc wireless sensor networks, a large number of applications require location awareness of the sensors, including tracking, environmental monitoring and many military applications. Without the knowledge of its position, the information captured by a sensor becomes obsolete. The main building block of these networks is a low-cost sensor, with low power resources, leaving no room to (absolute) selfpositioning device. To overcome this drawback, one includes in the network a small number of sensors with known positions (and sometimes high power and communication capabilities). These sensors, often known by anchors or beacons and designated hereafter by the latter, communicate to other sensors information allowing the latter to estimate their positions. For this purpose, each sensor determines ranging (distance) measurements with other sensors, from some measurements such as the received signal strength indication (RSSI), the connectivity, the hop count, the time difference of arrival, ... Most work on localization in sensor networks considers either multidimensional scaling (MDS) techniques or semidefinite programming (see [START_REF] Bachrach | Localization in sensor networks[END_REF][START_REF] Mao | Wireless sensor network localization techniques[END_REF], and references therein), in order to determine a function that links the ranging of the sensors to their positions, based on the known positions of some beacons.

{FirstName.SecondName}@utt.fr Introduced by Aronszajn in the mid 50s, the usefulness of reproducing kernels has been demonstrated in the last 15 years in the field of pattern recognition with the statistical learning theory and the so called kernel machines, such as support vector machines (SVM) and kernel principal component analysis (kernel-PCA) [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF]. Reproducing kernels provide new insights into sensor networks research field. This has been known for a while, as many researchers in sensor networks focus on detection, tracking and classification using kernel machines. In recent years, there has been an increasing interest in this framework for localizing sensors, with kernel-PCA [START_REF] Jenkins | Relative localization from pairwise distance relationships using kernel PCA[END_REF][START_REF] Essoloh | Distributed localization in wireless sensor networks as a preimage problem in a reproducing kernel hilbert space[END_REF], SVM [START_REF] Nguyen | Learning in decentralized systems: A nonparametric approach[END_REF], least squares regression [START_REF] Zhu | Ad hoc sensor network localization using distributed kernel regression algorithms[END_REF], and manifold regularization [START_REF] Pan | Adaptive localization in a dynamic wifi environment through multi-view learning[END_REF].

In this paper, we derive a two-stage strategy. First we seek a mapping function between the ranging and the inner products between positions of a given sensor with beacons. Learnt with data available from beacons, it is then applied to any sensor, leading to an estimate of the inner products between its (unknown) position and the (known) positions of the beacons. In the second stage, we determine the position of the sensor from these estimated inner products. It turns out that the first stage is a matrix completion problem, where the inner-product matrix is completed from the (entirely available) ranging matrix, and thus can be solved with the recently introduced matrix regression method [START_REF] Yamanishi | Kernel matrix regression[END_REF]. By learning the regression from available data in both matrices simultaneously, thus from inter-beacon information, we show that this reduces to a linear optimization problem. The second stage can be solved by considering the Nyström method, a technique for approximating kernel matrices in the machine leaning community. We investigate a distributed version of the method, by solving locally the optimization problem. We emphasize that the proposed method is independent of the ranging type, thus can be applied to RSSI, hop count, or any other ranging information. This paper is organized as follows. We begin by presenting the matrix regression method of sensor localization. Then, in section 3, we derive a technique for determining the position of the sensors. We propose, in section 4, a distributed algorithm, taking into account only nearby beacons for the considered sensor. Finally, computer simulations are carried out to validate the proposed approach.

THE MATRIX REGRESSION METHOD

Consider a network of m sensors of unknown positions and n beacons of known positions, living in a d-dimension (2D or 3D) space. Let X and Y be the coordinate matrices of beacons and sensors, respectively, of size n-byd and m-by-d, and [X ⊤ Y ⊤ ] ⊤ the overall coordinate matrix. The inner product between their coordinates is given by P

= X Y [X ⊤ Y ⊤ ],
which can be decomposed into four block submatrices, P x = XX ⊤ , P yx = YX ⊤ , P xy = P ⊤ yx , and P y = YY ⊤ , as illustrated in [START_REF] Bachrach | Localization in sensor networks[END_REF] with unknown submatrices set to gray-color. On the other hand, we have the overall ranging matrix, denoted by K with entries κ(x i , x j ), similarly decomposed into K x , K yx , K xy , and K y , as given in [START_REF] Bachrach | Localization in sensor networks[END_REF].

K x K xy K yx K y K → P x P xy P yx P y P (1)
In a conventional regression problem, one seeks a function φ(•) that links an input variable x into a response (output) variable z, under the constraints φ(x i ) = z i for all available training data {(x 1 , z 1 ), . . . , (x n , z n )}. While there exists an infinite number of functions verifying such constraints, one considers functions with some regularizing properties (such as smoothness). This can be done by restricting the hypothesis space to the RKHS of a given reproducing kernel, say κ(•, •). Moreover, from the Representer Theorem [START_REF] Schölkopf | A generalized representer theorem[END_REF], the optimal function has the form

φ(•) = n k=1 α k κ(x k , •). (2)
For instance, this is true for kernel-PCA, where each principal axis φ(•) is determined by its n coefficients α k , obtained by an eigen-decomposition of the n-by-n matrix of entities κ(x i , x j ), thus K x . Since φ(x) corresponds to the principal coordinate of x, the latter can be represented into a lowdimensional space by considering only a couple of principal coordinates. Since this is the essence of both MDS and kernel-PCA techniques, localization in sensor networks using kernel-PCA is proposed in [START_REF] Jenkins | Relative localization from pairwise distance relationships using kernel PCA[END_REF], or more recently [START_REF] Essoloh | Distributed localization in wireless sensor networks as a preimage problem in a reproducing kernel hilbert space[END_REF] (see [START_REF] Williams | On a connection between kernel PCA and metric multidimensional scaling[END_REF] for a connection to MDS).

In what follows, we consider the general case of determining a set of optimal functions, fully described by their coefficients, which identifies the mapping between the two matrices described in [START_REF] Bachrach | Localization in sensor networks[END_REF]. This is known as a matrix regression problem [START_REF] Yamanishi | Kernel matrix regression[END_REF], between the input data κ(x i , x j ) of matrix K and the output x i x ⊤ j of P. We learn this problem from the available input-output couples, i.e. K x and P x . For this, we consider a model of the form

Ψ(x i , x j ) = x i x ⊤ j , (3)
and determine it from inter-beacon information. As above, we consider a particular form of Ψ, by representing each x in some coordinates obtained with a set of functions,

[φ 1 (x) φ 2 (x) • • • ] ⊤ , leading to Ψ(x i , x j ) = [φ 1 (x i ) φ 2 (x i ) • • • ][φ 1 (x j ) φ 2 (x j ) • • • ] ⊤ = h φ h (x i )φ h (x j ).
In analogy to kernel-PCA where these are principal coordinates obtained from principal axes constructed from available data κ(x k , x ℓ ), we consider the same form as (2) for all φ h 's, constructed from K x . Therefore one should determine for each φ h the optimal coefficient vector

α = [α 1 α 2 • • • α n ] ⊤ , with φ h (x i ) φ h (x j ) = n k=1 α k κ(x k , x i ) n ℓ=1 α ℓ κ(x ℓ , x j ) = κ ⊤ i α α ⊤ κ j
where κ i is the i-th column vector of K x . From the sum of these terms and (3), we get x i x ⊤ j = Ψ(x i , x j ) = κ ⊤ i A κ j , where A is a coefficient matrix. Since this should be satisfied for all beacons, i.e. i, j = 1, . . . , n, we can write the optimization problem in matrix form, with

min A P x -K ⊤ x A K x 2 F , (4)
where • F is Frobenius norm. Once we obtain the optimal coefficient matrix A, we can apply the resulting map to sensors with unknown positions, with x i y ⊤ j = Ψ(x i , y j ). From expressions above, we obtain the matrix expression

P xy = K ⊤ x A K xy . (5)
The optimization problem (4)-( 5) can be solved by writing (4) as a generalized eigen-decomposition problem and injecting the resulting matrix in [START_REF] Essoloh | Distributed localization in wireless sensor networks as a preimage problem in a reproducing kernel hilbert space[END_REF]. As we notice that both expressions contain the matrix T = K ⊤

x A, we propose to solve the following equivalent optimization problem :

min T P x -T K x 2 F , and P xy = T K xy . (6)
This is a linear problem yielding T = P x K -1

x , and thus

P xy = P x K -1 x K xy . (7)

SENSOR POSITION ESTIMATION

After estimating the matrix P xy of inner products of positions between beacons and sensors, one has to find the coordinates of the latter. For this purpose, we take advantage of the Nyström method, initially developed in the machine learning community to approximate a matrix by another matrix of lower rank [START_REF] Williams | Using the Nyström method to speed up kernel machines[END_REF]. In our case, on the one hand we have by construction P x = XX ⊤ , thus a d-rank matrix. From its eigen-decomposition, we have

P x = U d Λ d U ⊤ d ,
where Λ d is a diagonal matrix of the d nonzero eigenvalues of P x , and U d the matrix whose columns are the corresponding eigenvectors. By identification, we get

X = U d (Λ d ) 1/2 (8)
On the other hand, we can write P xy = XY ⊤ where Y is the coordinate matrix to be identified. By injecting [START_REF] Pan | Adaptive localization in a dynamic wifi environment through multi-view learning[END_REF] into this definition, we get the coordinates of the m sensors from

Y ⊤ = (Λ d ) -1/2 U ⊤ d P xy . ( 9
)
Since the resulting coordinates are determined in the space defined by the eigenvectors, one must conduct a final step of mapping them, with a linear (or affine to be more precise) transformation, into the initial space of beacons. Such step is commun to MDS techniques.

DISTRIBUTED ALGORITHM

In section 2, we considered completing the inner-product matrix by inverting the ranging matrix of the beacons. This stipulates that the beacons communicate to each other, in a peer-to-peer fashion or a multi-hop strategy. However, in practice, beacons may not be in the range of each other. Moreover, the matrices may become too cumbersome to invert and manipulate for a large scale sensor network. To overcome this problem, we propose a distributed algorithm, where each sensor gets information from nearby beacons in order to find its own position. In other words, any sensor i defines a set of neighboring 1 beacons with a submatrix of K x and its counterpart in P x , denoted respectively K i and P i , and X i the corresponding coordinates. While the global optimization problem (6) leads to expression [START_REF] Zhu | Ad hoc sensor network localization using distributed kernel regression algorithms[END_REF], by considering the distributed approach, we get

p i = P i K -1 i κ i , (10)
where p i is the inner-product column vector of positions between sensor i and its nearby beacons, and κ i the column vector of ranging between them.

The corresponding coordinates of this sensor can be revealed by writing locally the expression (9), obtained from 1 Different strategies can be proposed to define the neighborhood of a sensor. This can be done by examining the ranging values, where high values correspond to neighbors. We fix their number in simulations.

for each sensor i Find the nearby beacons [dump,ind]=sort(K(i,1:n)) Consider the closes nc beacons ind=ind(end:-1:end-nc+1) Get ranging between these beacons Ki=K(ind,ind) Get inner products between them Pi=Xn(ind,:) * Xn(ind,:)' Consider ranging with them ki=K(i,ind) Compute inner products with [START_REF] Schölkopf | A generalized representer theorem[END_REF] pi=Pi * inv(Ki) * ki Determine position with [START_REF] Williams | On a connection between kernel PCA and metric multidimensional scaling[END_REF] y=pi/(Xn(ind,:)')

Table 1. Pseudocode of the distributed algorithm.

an eigen-decomposition problem. Then, a mapping transformation must be carried out as presented above, by considering this time only neighboring beacons. While this becomes fairly cumbersome for each sensor, we propose an alternative approach to find the coordinates, based on the pseudoinverse operator. For this purpose, we rewrite the problem as the following optimization problem min

y p i -X i y ⊤ 2 F .
It is well known that the solution of this linear system is given by the left pseudoinverse of the matrix X i , with

y = X ⊤ i X i -1 X ⊤ i p i . (11)
We emphasize on the fact that we don't need to apply a mapping to localize the sensors with respect to beacons. The simplicity of the algorithm is illustrated in Table 1.

SIMULATIONS

To illustrate our method, we consider a configuration similar to the one proposed in [START_REF] Patwari | Manifold learning algorithms for localization in wireless sensor networks[END_REF], with ranging between two sensors is only a function of the distance between them, with κ(x i , x j ) = expx ix j 2 /2σ 2 , where σ is a parameter corresponding to the range of the sensors. Next, sensors are randomly spread on a 1-by-1 square region. In a first series of experiments, we apply the centralized algorithm to a network of 200 sensors, and study the influence of both the number of beacons and the range parameter. For this, we take n = 3, 4, . . . , 20 and σ = 0.1, 0.2, . . . , 0.9. In Fig. 1, we plot the resulting rootmean-square error, averaged over 100 trials. As expected, the localization error decreases as the number of beacons increases, and the visibility between sensors is high. By taking for instance 15 beacons (almost 7% of the sensors), and σ = 0.75, we get the topology illustrated in Fig. 2.

In a second experiment, we consider a large scale network of 1000 sensors of low range, with σ = 0.3. We deploy also 50 beacons with the same characteristics. This setting results in inverting and manipulating large sparse matrices, since there is low visibility between these entities. For these reasons, we consider the distributed algorithm, each sensor determines its coordinates with information from the 5 closest beacons. Fig. 3 illustrates the resulting topology, with a root-mean-square error of 0.018. Fig. 3. Topology constructed by the distributed algorithm for a large scale network (same legend as Fig. 2).

4 NFig. 1 .

 41 Fig. 1. Root-mean-square error on positions with the centralized algorithm, as a function of n and σ.

Fig. 2 .

 2 Fig. 2. Topology constructed by the centralized algorithm. The beacons are represented by , sensors positions by +, and their estimations by •.

CONCLUSION

In this paper, we took advantage of recent works in kernel machines for solving the localization problem in sensor networks. We showed that the matrix regression method allows us to estimate unknown positions of sensors. We derived a distributed algorithm, based on information from local neighborhood of each sensor. There are several directions for further research, including mobile ad hoc network (MANet), with an iterative update of the coordinates.