Paul Honeine

Farah Mourad

Maya Kallas

Hichem Snoussi

Hassan Amoud

Clovis Francis

WIRELESS SENSOR NETWORKS IN BIOMEDICAL: BODY AREA NETWORKS

The rapid growth in biomedical sensors, low-power circuits and wireless communications has enabled a new generation of wireless sensor networks: the body area networks. These networks are composed of tiny, cheap and low-power biomedical nodes, mainly dedicated for healthcare monitoring applications. The objective of these applications is to ensure a continuous monitoring of vital parameters of patients, while giving them the freedom of motion and thereby better quality of healthcare. This paper shows a comparison of body area networks to the wireless sensor networks. In particular, it shows how body area networks borrow and enhance ideas from wireless sensor networks. A study of energy consumption and heat absorption problems is developed for illustration.

INTRODUCTION

Nowadays, biomedical sensors are used solely in health monitoring applications, as with the electrocardiograms. In other words, they are not yet integrated into communicating networks. However, with the increase of the vital body signs to be monitored simultaneously, many biomedical sensors need to work together to collect measurements at the same time. Moreover, it is often important to communicate the measured data from a sensing device to others. Thus, it becomes crucial to have a network of communicating biomedical sensors. This would have a radical impact on the quality of life of patients and on their treatment success rates. It could have a wide range of future applications, as well, such as monitoring of cardiovascular diseases, diabetes and asthma, consulting cases via telemedicine or health-care systems, etc.

The Body Area Networks (BANs) have recently emerged as a by-product to do so. They are a new generation of Wireless Sensor Networks (WSNs), suited to monitor human bodies. A BAN is composed of a collection of tiny nodes, equipped with biomedical sensors, motion detectors and wireless communication devices. These nodes collect vital body This work was partly supported by the Franco-Lebanese CEDRE program.

signs and then transmit them in a wireless manner to a central unit where all collected information is processed. Due to their wireless nature, BAN nodes have many advantages, such as omnipresence connectivity, mobility and inter-operability [START_REF] Cypher | Prevailing over wires in healthcare environments: benefits and challenges[END_REF]. Some nodes may also be equipped with actuators, such as pacemakers or instruments to store and inject medicines.

In this paper, we review some issues in BANs. We show how BAN borrows and enhances ideas from the WSN literature. More specifically, we study the energy consumption and the heat absorption problems, both problems being tightly related.

FROM WSN TO BAN

Wireless sensor networks

Recent advances in miniaturization and wireless communication have enabled a new class of tiny powerful computers: Wireless Sensor Networks (WSNs). WSNs consist of largescale, highly distributed and self-organized networks, mainly developed for surveillance and monitoring. They are composed of a large number of tiny and cheap nodes, having computing, sensing and communication capabilities. WSNs were initially meant for military applications such as large battlefield surveillance. In these applications, WSNs were assumed to have cubic-millimeter computers, called smart dust [START_REF] Warneke | Smart dust: Communicating with a cubic-millimeter computer[END_REF]. To the best of our knowledge, the smart dust concept has not been realized. However, the issues and challenges of WSNs pave the way to many new opportunities. It is thus expected to use WSNs in a broad range of applications, including environmental monitoring, animal migration monitoring, security surveillance, industrial control, as well as personal healthcare monitoring, etc.

Description of body area networks

The technology behind the WSNs is still under development. However, the research on WSNs does not deal with the challenges associated with human body 1 A BAN is a network of wearable on-body computing devices. It can also includes in vivo implanted biosensor devices. According to the IEEE 802.15 (Task Group 6), a BAN consists of "low power devices operating on, in or around the human body (but not limited to humans) to serve a variety of applications including medical, consumer electronics / personal entertainment and other". BANs have a wide range of applications, as, for instance, in gaming or virtual reality. However, the main challenges in terms of research or engineering remain in the biomedical and healthcare monitoring applications. Indeed, evolution of BANs should follow the ever-increasing development in the medical domain; its main objective being to ensure constant vigilant and pervasive monitoring of patients at home or at work.

Requirements in body area networks

Although BANs are essentially WSNs, they take the WSNs to their extreme in many directions. Indeed, BANs for instance are limited to the human body where very few nodes are deployed. Having few sensing devices, data loss in BANs could be significant, as opposed to WSNs, where nodes yield redundant information. Therefore, each node needs to provide accurate measurements, especially for nodes giving vital information, such as electrocardiogram (ECG) measurements. This may require additional measurements to ensure quality-of-service and real-time data delivery, since reliability of measurements is a must in the medical domain. Moreover, BANs often consist of networks with devices performing diverse tasks. Examples of these devices include ECG, blood pressure, pulse oximetry, temperature, respiratory, etc. Such heterogeneous characteristics require different frequency acquisition rates as well as different transmission rates.

On the other hand, some devices of the BAN need to be implanted, in vivo, inside the human body. Implanted biosensors should be biocompatible, robust, as replacement is diffi-cult, and with low energy requirements, since energy is hard to supply. In addition, in BAN, nodes are in close proximity to or inside the human body. This may lead to high absorption of the electromagnetic radio-frequency waves, and thus to an increase of the body temperature. According to the safety limits for exposure to radio-frequency energy, communications should be uniformly distributed among all the nodes. Moreover, as the human body is not rigid, frequent changes in nodes positions should be considered in the network topology. For these reasons, a special attention should be paid to the design of the BAN, including the routing strategy.

ENERGY SCAVENGING

The energy consumption in BANs is crucial, especially in implanted biosensors, since they are inaccessible and difficult to replace. For instance, a pacemaker battery usually lasts for 5 to 10 years. Since the battery is sealed inside the pacemaker, replacing the battery leads to the replacement of the entire system. In order to increase the lifetime of battery-powered devices, extra amount of energy can be delivered by energy harvesting, i.e., energy scavenging [START_REF] Paradiso | Energy scavenging for mobile and wireless electronics[END_REF]. As opposed to WSNs, where a great attention is paid to solar energy, BANs introduced a new nature of energy scavenging. Devices can scavenge power from the human body using body heat and body vibrations. Although heat and motion of the human body bring new opportunities, they also introduce new challenges, as shown in the following.

Body motion

Electrical power can be generated from motion and vibration, with the so-called inertial-power scavenging. Natural human body motion may be easily converted into electrical power. This is illustrated by electronic self-winding watches based on moving mechanics, such as the ETA Autoquartz and the Seiko Automatic Generating System. More promising approaches are based on piezoelectric or capacitive generators. An example of electric power generators is the electric shoes [START_REF] Antaki | A gait-powered autologous battery charging system for artificial organs[END_REF][START_REF] Shenck | Energy scavenging with shoemounted piezoelectrics[END_REF]. For instance, authors of [START_REF] Antaki | A gait-powered autologous battery charging system for artificial organs[END_REF] propose electric shoes to provide autonomy to artificial organs.

While motion may provide power to BANs, it introduces new challenges. In fact, as the human body is not rigid, thus the devices can move, at least relatively to each others. Thus, frequent changes in the network topology should be considered in the design of the network architecture. BANs should be robust to these changes, for instance by adapting the routing strategy, and the self-organization of the nodes.

Body heat

Temperature difference can be directly converted into electrical power, as given by the Peltier-Seebeck effect. As any heat engine, the efficiency of thermoelectric energy is limited by the Carnot's theorem. In practice, large thermal gradients are required to create satisfactory levels of voltage and power. This limits the application of thermoelectric energy in BANs, provided by body heat over ambient temperature. Nevertheless, recent works are very encouraging [START_REF] Stevens | Optimal design of small ∆T thermoelectric generation systems[END_REF], as illustrated by the Seiko Thermic watch.

While environment temperature cannot be controlled, human body temperature is essentially 3 37 • C (98.6 • F). In BANs, the temperature is likely to be different than this common value. In fact, radio-frequency transmission of wireless nodes located on-body or implanted contribute to radiation absorption, which may result into thermal effects. For safety reasons, the radiation absorption should be reduced, by controlling the Specific Absorption Rate (SAR) [START_REF] Ren | Rate control to reduce bioeffects in wireless biomedical sensor networks[END_REF] and using an appropriate routing protocol (see Section 4.2).

ROUTING INFORMATION: FROM WSN TO BAN

In this section, the routing problem in BANs is tackled with techniques from the WSNs literature. Our aim is to show how BANs borrow and enhance ideas from WSNs. Three directions are presented here, with essentially the same goal: reducing the absorption of power in the human body, which is dissipated as heat. Such a problem is addressed by a clusterbased routing, a temperature-aware routing, and using an onbody node for implanted networks.

Routing by clustering

In WSNs, several properties for routing strategies are desirable: maximize system lifetime by uniform battery usage of nodes, and application to large-scale networks. Such properties are not satisfied by classical network protocols, including the direct transmission, multi-hop, and clustering.

The LEACH (Low Energy Adaptive Clustering Hierarchy) [START_REF] Heinzelman | Energyefficient communication protocol for wireless microsensor networks[END_REF] strategy fulfills the desired properties by combining, on the one hand a hierarchical protocol with clusterhead nodes to aggregate and compress the data, and on the other hand a randomized rotation of the clusterhead nodes. In fact, a cluster-based hierarchy allows small communication distances for most nodes, the non-clusterhead nodes [START_REF] Younis | Node clustering in wireless sensor networks: Recent developments and deployment challenges[END_REF]. As opposed to conventional clustering strategies, LEACH uses 3 Body temperature depends upon the place at which the measurement is made (slightly higher when taken internally), and fluctuates over the day. rotating election of clusterhead nodes and adapting the corresponding clusters, which allows the energy requirements to be uniformly distributed among all the senors. Moreover, data is compressed at each clusterhead node, accomplishing an additional reduction in the energy dissipation.

Several algorithms were derived for BAN based on LEACH, using a node clustering structure. For instance, the authors of [START_REF] Watteyne | Anybody: a self-organization protocol for body area networks[END_REF] focus on the self-organization of nodes, as opposed to the randomness election in LEACH which may lead to some isolated nodes. In [START_REF] Culpepper | Design and analysis of hybrid indirect transmissions (hit) for data gathering in wireless micro sensor networks[END_REF], the authors combine, on the one hand clustering to reduce the number of long direct communications, and on the other hand parallel, multi-hop indirect transmissions, the co-called chains. Such combination has significantly better energy efficiency than LEACH, as illustrated in [START_REF] Moh | On data gathering protocols for in-body biomedical sensor networks[END_REF] for in-body biomedical sensor networks. It is worth noting that LEACH is optimized for large-scale networks, which is not the case for BAN.

Temperature-aware routing

The Specific Absorption Rate (SAR) is defined as power absorbed per unit mass of the tissue, a measure in determining the amount of power lost due to heat dissipation. In BAN, the SAR should be taken into account, either locally or globally using its average over the entire mass of tissue between the transmitter and the receiver. Because of the SAR consideration, the thermal effects of BAN should be considered in the routing protocol. This can be done by balancing the communication over the nodes, as illustrated above with the LEACH algorithm. Moreover, one may also consider the temperature to select appropriate routing paths.

Temperature-aware routing was initially investigated in [START_REF] Tang | Communication scheduling to minimize thermal effects of implanted biosensor networks in homogeneous tissue[END_REF], by rotating clusterhead nodes based on leadership history and node location. In [START_REF] Bag | Energy efficient thermal aware routing algorithms for embedded biomedical sensor networks[END_REF], an (adaptive) least temperature routing is proposed by routing data away from high temperature spots and area. To this end, multi-hop routing is applied, where each node selects its neighboring node with the lowest temperature as its next hop. By selecting the next hop, such approach is not optimized globally. To overcome this drawback, one may consider the least total temperature route, by including the shortest path routing. In [START_REF] Takahashi | Temperature-aware routing for telemedicine applications in embedded biomedical sensor networks[END_REF], node temperature is converted into graph weights and minimum temperature route is obtained using Dijkstra's algorithm.

Routing: in vivo implanted networks

More constraints are imposed by in vivo implanted sensor networks, including restrictions due to absorption of communication power in the tissue (SAR), which is dissipated as heat. For this reason, classical air-as-medium propagation models, where losses are due to freespace wave propagation and multipath fading, are no longer appropriate. In the human body, the saline-water nature of the human tissues contributes considerably to the absorption of the electromagnetic waves in communication. There has been considerable research going on in the field of measurement of calculation of SAR in human tissues, and specification of safe absorption rates [START_REF]A practical guide to the determination of human exposure to radiofrequency fields[END_REF].

It is well known that the in-body propagation loss is much more than that of on-body propagation. Thus, routing data among implanted nodes may be inappropriate for long range communications. In [START_REF] Zhen | A body surface coordinator for implanted biosensor networks[END_REF], a new network architecture is proposed for the implanted networks by introducing an on-body coordinator. The coordinator plays the role of a beacon, by forwarding the data from one implanted node to another over long distance more safely and efficiently. Since a coordinator is on-body, it may include other functionalities, such as monitoring and more power resources.

FINAL REMARKS

In this paper, our aim was to show how BAN and WSN are intimately connected, borrowing and enhancing ideas from each others. This was illustrated on the routing problem, where BAN incorporates additional constraints including energy consumption, SAR and temperature considerations. All these constraints are intimately connected. It is obvious that this paper is not a survey on BAN. For instance, the physiological biosensors and the MAC layer were not studied, nor security issues [START_REF] Stuart | Privacy and security in biomedical applications of wireless sensor networks[END_REF][START_REF] Hanlen | Key-sharing via channel randomness in narrowband body area networks: is everyday movement sufficient?[END_REF]. We refer the interested reader to more comprehensive surveys [START_REF] Ullah | A comprehensive survey of wireless body area networks[END_REF][START_REF] Chen | Body area networks: A survey[END_REF][START_REF] Latré | A survey on wireless body area networks[END_REF]. It is worth noting that most of these surveys were published in 2011, subsequent to the initial submission of this work.

Table 1 .

 1 Some differences between WSNs and BANs.

	monitoring. For this rea-

Body area networks are not limited to humans.

BANs are also known as wireless body area networks, or body sensor networks.


Zineb Noumir

Paul Honeine

Cédric Richard

MULTI-CLASS LEAST SQUARES CLASSIFICATION AT BINARY-CLASSIFICATION COMPLEXITY

This paper deals with multi-class classification problems. Many methods extend binary classifiers to operate a multiclass task, with strategies such as the one-vs-one and the onevs-all schemes. However, the computational cost of such techniques is highly dependent on the number of available classes. We present a method for multi-class classification, with a computational complexity essentially independent of the number of classes. To this end, we exploit recent developments in multifunctional optimization in machine learning. We show that in the proposed algorithm, labels only appear in terms of inner products, in the same way as input data emerge as inner products in kernel machines via the so-called the kernel trick. Experimental results on real data show that the proposed method reduces efficiently the computational time of the classification task without sacrificing its generalization ability.

INTRODUCTION

Many challenging classification tasks are multi-class problems, often shown to be more difficult than the binary classification. They imply classifying a data observation into one of a set of possible classes. Multi-class problems are encountered in most applications in signal and image processing, for instance for optical character recognition, speech application [START_REF] Klautau | Combined binary classifiers with applications to speech recognition[END_REF], face recognition [START_REF] Lihong | Face recognition based on multi-class SVM[END_REF], and classification of urban structures in an image [START_REF] Honeine | The angular kernel in machine learning for hyperspectral data classification[END_REF].

Multi-class classification methods can be roughly divided into two categories. The first one consists of a single machine strategy, where a single optimization problem is solved in order to determine the multi-class classifier. These techniques require specific optimization algorithms, often with high computational cost [START_REF] Weston | Support Vector Machines for Multi-Class Pattern Recognition[END_REF][START_REF] Crammer | On the algorithmic implementation of multiclass kernel-based vector machines[END_REF]. The second category attempts to take advantage of the performance of binary classifiers such as state-of-the-art least squares classifiers and support vector machines. To this end, the multi-class classification problem is decomposed into a number of binary classification subproblems. The goal is to solve these subproblems, then combine This work was supported by ANR-08-SECU-013-02 VigiRes'Eau.

their results to determine the multi-class solution (see for instance [START_REF] Suykens | Multiclass least squares support vector machines[END_REF]). Two strategies exist to extend binary classifiers to multi-class, the one-versus-all strategy where binary classifiers are constructed to separate one class from all the other classes, and the one-versus-one when separating one class from another (see [START_REF] Rifkin | In defense of one-vs-all classification[END_REF] and references therein). For a problem of m classes, these strategies require solving and combining m binary classifiers for the former, and m(m -1)/2 for the latter. From the quantitative and comparative study given in [START_REF] Rifkin | In defense of one-vs-all classification[END_REF] with least squares binary classifiers, the one-vs-all strategy gives at least similar, if not better, results than other methods, including the single machines. Still, this requires m binary classifiers, each constructed using the whole available training data.

In this paper, we show that one can solve a multi-class classification problem, with essentially the same computational complexity as a binary classifier. To this end, we take advantage of recent work in multifunctional optimization in machine learning [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF]. We recast the multi-class problem as the estimation of a vector of outputs defining the class membership. By using a least squares formalism, the resulting optimization problem is roughly similar to the one given by a binary classifier. It turns out that its computational complexity is essentially independent of the number of classes, in the same sense as the computational cost in kernel-machines is independent of the dimensionality of the input data, thanks to the so-called kernel trick.

The rest of the paper is organized as follows. Section 2 outlines the least squares approach for a binary classification problem. We describe the proposed multi-class least squares algorithm in Section 3. Section 4 illustrates results obtained with our algorithm, with an image classification problem. Conclusions and further directions are given in Section 5.

LEAST SQUARES BINARY CLASSIFICATION

In supervised learning, one seeks a function that predicts well some output from a given input, based on a set of training input-output data, (x k , y k ), for k = 1, 2, . . . , N . To measure the excellence of the learned function f (•), a loss function is considered, such as the commonly used square error, namely,

ℓ(f (x k ), y k ) = |f (x k ) -y k | 2 .
The minimization of this square error is considered for binary classification problems in [START_REF] Rifkin | Everything old is new again: A fresh look at historical approaches in machines learning[END_REF][START_REF] Zhang | SVM vs regularized least squares classification[END_REF] (with connections to the Fisher discriminant analysis, however beyond the scope of this work). In this case, the output corresponds to the label, in practice y k = ±1.

While this loss function is minimized over all the training data, one also imposes regularity to the function, using for instance the Tikhonov regularization [START_REF] Tikhonov | Solutions of Ill-Posed Problems[END_REF]. This gives the following least squares optimization problem:

min f N j=1 |f (x j ) -y j | 2 + γ f 2 ,
where γ is some positive tradeoff parameter. Moreover, in a kernel-based formalism [START_REF] Schölkopf | A generalized representer theorem[END_REF], one seeks a function of the form

f (•) = N j=1 α j y j κ(x j , •), (1)
with κ(•, •) being a positive semi-definite function, such as the linear kernel κ(x i , x j ) = x T i x j and the Gaussian kernel κ(x i ,

x j ) = exp( 1 2σ 2 x i -x j 2
) with some bandwidth parameter σ. It is work noting that, since one often sets y j = ±1 for a binary classification task, one may inject its sign into α j with β j = y j α j , and take the equivalent expansion f (•) = N j=1 β j κ(x j , •) (see for instance [START_REF] Rifkin | Everything old is new again: A fresh look at historical approaches in machines learning[END_REF]). By substituting this expansion into the least squares optimization problem, we obtain a classical matrix formulation, which can be solved using off-the-shelf linear algebra techniques. The resulting coefficients are used in the decision for any x, by comparing f (x) to the threshold, often set to zero, according to the rule

f (x) = N j=1 α j y j κ(x j , x) y=+1 ≷ y=-1 0.
Binary classifiers can be easily extended to operate for a multi-class classification task. Let m be the number of classes. In a one-vs-one strategy, one constructs m(m -1)/2 binary classifiers by taking only data from each pair of classes, while in a one-vs-all strategy, m binary classifiers are considered with all training data. Each binary classifier is defined by the decision function, of the form

f k (•) = N j=1 α j,k y j κ(x j , •). (2)
Therefore, one needs to estimate m × N coefficients for the one-vs-all strategy. Next, we propose to estimate only N unknown variables by imposing some relation between these functions, f 1 (•), f 2 (•), . . . , f m (•).

MULTI-CLASS LEAST SQUARES CLASSIFICATION

In a multi-class classification problem, we consider a set of N training data, belonging to any of the m available classes. We propose to encode1 the class membership using a m-column label vector y k , with its t-th entry given by [y k ] t = 1 if data x k belongs to class t; 0 otherwise for t = 1, 2, . . . , m. We propose to estimate all the functions in (2), using a single optimization problem, by requiring some relation between them. We regroup the m functions in the form

f (•) = [f 1 (•) f 2 (•) • • • f m (•)] T .
By analogy with (1) and following [START_REF] Szedmak | Multiclass learning at one-class complexity[END_REF], we propose the following expansion for f (•)

f (•) = N j=1
α j y j κ(x j , •).

In the proposed expression, all the functions in f (•) share the same scalar value α j κ(x j , x), for some j. This allows us to have only N unknowns, the α j 's, in the same way as a single binary classifier. Next, we show how to estimate these coefficients.

In a least squares sense, we consider the following optimization problem

min f N j=1 f(x j ) -y j 2 + γ R(f ),
where the regularization term is given by

R(f ) = N i,j=1
α i α j y T i y j κ(x i , x j ).

By substituting the expression of f (•) in this optimization problem, we get

min α N i=1 y T i y i -2 N i=1 y T i N j=1 α j κ(x i , x j ) y j + N i=1 N j,k=1 α j α k κ(x i , x j )κ(x i , x k ) y T j y k + γ N i,j=1 α i α j y T i y j κ(x i , x j ).
We drop the first term, since it is independent of α. In matrix form, this optimization problem can be written as

min α -2 d α + α T Gα + γ α T Ωα, (3)
where G is a matrix whose (j, k)-th entry is

[G] j,k = y T j y k K 2 (j, k)
d is a vector whose j-th entry is

[d] j = N i=1 κ(x i , x j ) y T i y j ,
and Ω is a matrix whose (i, j)-th entry is

[Ω] i,j = y T i y j κ(x i , x j ).

By taking the gradient of the objective function in (3) with respect to α, namely

-d + Gα + γ Ω α,
and setting it to zero, we obtain the final solution

(G + γ Ω) α = d. ( 4
)
We conclude that a multi-class regularized least squares classifier can be obtained by solving a single system of N linear equations with N unknowns. This requires the inversion of a N -by-N matrix, G + γ Ω. The computational complexity of this algorithm is cubic in the number of training data, N , but independent of the number of classes m. This is made possible here thanks to the fact that the label vectors only appear in terms of inner products, with y T i y j in Ω, G and d. This reflects an analogy with the kernel trick in kernel machines, where data are involved only in terms of inner products, namely x T i x j . To predict the class membership of any x, the decision rule compares the output f (x) with the set of label vectors y 1 , y 2 , . . . y m , with

d(x) = arg max t y T t f (x).
This decision rule can also be written using inner products between label vectors, namely

d(x) = arg max t N i=1 α i y T t y i κ(x i , x). (5)

EXPERIMENTATIONS

Two experiments were conducted to illustrate the pertinence of the proposed approach. First, an illustrative application is considered, with the IRIS data from the UCI Machine Learning Repository, often used as a benchmark for classification algorithms. It consists of samples of flowers representing m = 3 iris species. Each species consists of 50 observations, four features were measured from each sample, they are the length and the width of sepale and petale. For illustration, only sepale length and petale width were used. A set of N = 120 training data was used for learning the multi-class classifier. Some preliminary experiments were conducted to tune the parameters, leading to γ = 100 and σ = 0.60 for the Gaussian kernel. Figure 1 illustrates the resulting separating boundaries, and shows the decision rule given in [START_REF] Crammer | On the algorithmic implementation of multiclass kernel-based vector machines[END_REF] in the third dimension. For a comparative study, we considered the least squares one-vs-all strategy, as given in Section 2. Using the remaining set of 30 data, both multiclass classification methods gave essentially similar classification error, equal to 3%. In an attempt to provide a measure of computational requirements2 , the one-vs-all classifier was trained with a (average ) total CPU time of up to 0.011 seconds, while the proposed algorithm required only 0.003 seconds.

For the second application, we considered a multi-class classification task, where the number of available data is not 1.

the same for all classes. Recently adopted in remote sensing, hyperspectral images are cubes of data, measuring spectral composition within a spatial view. As opposed to the conventional 3-color system, the spectral information over a hundred of bands provides greater analysis of the composition of objects in an image scene. In monitoring urban structures with airborne or satellite images, one is often confronted with a large number of classes.

The hyperspectral image, provided by the HySenS project, is of the University of Pavia, Italy, with a 300-by-300 pixels and 103 frequency bands, illustrated in Figure 2 (left). Ground truth information about 7 classes were included to train and test the classifiers, as given in Table 1 illustrated in Figure 2 (middle and right) (see [START_REF] Honeine | The angular kernel in machine learning for hyperspectral data classification[END_REF] for more details). The Gaussian kernel was used, with its bandwidth value set to the maximum value in the training data. Both the proposed method and the one-vs-all method give almost similar results, with an overall error of 8% for both. It is worth noting that our algorithm estimates N = 1 492 coefficients by the inversion of a 1 492-by-1 492 matrix, independent of the number of classes, while the one-vs-all strategy estimates 7 × 1 492 coefficients for the m = 7 classes, thus operates the inversion of a 10 444-by-10 444 matrix.

CONCLUSIONS AND FUTURE WORK

In this paper, we studied the multi-class least squares classification, and showed that one can do multi-class classification at binary-classification complexity. While the idea is simple, its relevance was illustrated on real data, showing that our approach reduces efficiently the computational complexity without sacrificing its accuracy. In future work, we exploit further the proposed approach, by incorporating a nonlinear measure of similarity between label vectors, i.e., the use of a kernel on labels as opposed to the (linear) one in this paper.

Fig. 1 .

 1 Fig. 1. Illustration of the IRIS classification task, with the separating boundaries (-) and the decision rule taking values 1, 2 and 3.

Fig. 2 .

 2 Fig. 2. The hyperspectral image (slice at mid spectral-band) (left), with the spatial distribution of the training (middle) and test (right) data. The legend of the 7 classes is given in Table1.

Table 1 .

 1 The 7 classes with the ratio of train/test samples in the hyperspectral image.

	Class-name	#train	#test
	Asphalt	210	3 088
	Meadows	236	5 216
	Trees	196	1 279
	Metal sheets (painted)	256	1 345
	Bare soil	332	1 264
	Self-blocking bricks	225	1 693
	Shadow	28	215
		1 492 14 100

In this paper, we do not discuss the issue of optimal coding.

To offer a comparative study, both algorithms were implemented on a MATLAB running on a Windows.


Jie Chen

Jose Bermudez

Paul Honeine

A modified non-negative LMS algorithm and its stochastic behavior analysis

I. INTRODUCTION

Optimizing an objective function subject to constraints is very common in statistical signal and image processing. Constraints represent prior knowledge that may lead to proper solutions, or may reduce the feasible space of solutions when the problem is ill-conditioned. Non-negativity is probably one of the most commonly used constraints. It is often imposed on the parameters to estimate in order to avoid physically absurd and uninterpretable results. Non-negative least-square problems (NNLS) have been addressed in applications ranging from image deblurring [START_REF] Benvenuto | Nonnegative leastsquares image deblurring: improved gradient projection approaches[END_REF] to impulse response estimation [START_REF] Lin | Bayesian regularization and nonnegative deconvolution for room impulse response estimation[END_REF]. Non-negative matrix factorization (NMF) [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF], which is closely related to blind deconvolution problems, have found direct application in hyperspectral imaging [START_REF] Berry | Algorithms and applications for approximate nonnegative matrix factorization[END_REF].

A variety of methods have been proposed in the literature to tackle the NNLS problem. Active set algorithm of Lawson and Hanson [START_REF] Lawson | Solving Least Squares Problems[END_REF] is a batch resolution technique for NNLS problems, which has become a standard among the most frequently used methods. Projected gradient methods [START_REF] Calamai | Projected gradient methods for linearly constrained problems[END_REF], [START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF] are based on successive projections on the feasible region. Multiplicative algorithms are very popular to solve NMF problems [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF]. All these algorithms are however based on batch processing, which is not suitable for online system identification problems. In [START_REF] Chen | Non-negative least-mean-square algorithm[END_REF], [START_REF] Chen | System identification under non-negativity constraints[END_REF], online system identification subject to non-negativity constraints on the parameters to estimate was addressed. A LMS-type algorithm, called nonnegative LMS algorithm (NN-LMS) was proposed to compute solutions adaptively. It relies on stochastic gradient descent combined with a fixed point iteration scheme to converge to a solution that satisfies Karush-Kuhn-Tucker conditions.

During convergence of NN-LMS algorithm, it has been observed that weights may have different convergence rates and accuracy, especially those in the active set since they become smaller and finally tend to zero as iterations proceed. In this paper we introduce a modified NN-LMS algorithm in order to alleviate these unbalanced convergence rates. We also proposed analytical models to characterize the stochastic behavior of this algorithm.

II. MODIFIED NN-LMS ALGORITHM

Consider an unknown system, only characterized by a set of real-valued discrete-time responses to known stationary inputs.

The problem is to derive a transversal filter model

y(n) = α x(n) + z(n), (1)
with α = [α 1 , α 2 , . . . , α N ] the vector of model parameters, and

x(n) = [x(n), x(n -1), . . . , x(n -N + 1)
] the observed data vector. The input signal x(n) and the desired output signal y(n) are assumed zero-mean stationary. Sequence z(n) represents measurement noise and modeling errors. Due to the inherent physical characteristics of systems under investigation, the problem of identifying the optimum model can be formalized as follows

α o = arg min α J(α) subject to α i ≥ 0, ∀i, (2)
with

J(α) = E{[y(n) -α x(n)] 2 },
α o i [-∇ α J(α o )] i = 0, (3)
where the extra minus sign is just used to make a gradient descent of J(α) apparent. Applying a fixed-point iteration strategy and using usual stochastic gradient approximations, we obtain the NN-LMS algorithm [START_REF] Chen | Non-negative least-mean-square algorithm[END_REF], [START_REF] Chen | System identification under non-negativity constraints[END_REF]

α(n + 1) = α(n) + η e(n) D α (n) x(n) (4)
where η is a positive step size that must be tuned to construct a contraction scheme and to control the convergence rate. D α (n) stands for the diagonal matrix with diagonal entries given by α(n), and e(n

) = y(n) -α (n) x(n).
It is interesting to note that the α i (n)'s in the second term of the r.h.s of expression ( 4) are factors that modify the step size along each axis. This results in different convergence rates and accuracy for each component of α(n), especially for the weights corresponding to the active set as they become smaller through iterations and finally tend to zero. In order to compensate this unbalance of convergence rates, we now introduce the modified NN-LMS algorithm. Considering KKT condition (3), we can equivalently write

[α o i ] γ [-∇ α J(α o )] i = 0, (5)
with γ a positive real number. Implementing the fixed-point iteration strategy with equation ( 5) and using instantaneous estimates for the gradient leads us to the following algorithm

α(n + 1) = α(n) + η e(n) D x (n) α γ (n) (6)
Without ambiguity α γ (n) denotes the exponential value γ applied to each component of α(n). We recommend to specify parameter γ in the form γ = p q where p and q are both odd numbers, and 0 < p ≤ q. The oddness of p and q preserves the sign of α i (n) at each iteration. Although it can be shown that weights are always non-negative thanks to an appropriate step size, the fixed-point 0 can be reached by the left side 0 -with less strict step size selection. See [START_REF] Chen | Non-negative least-mean-square algorithm[END_REF] for more details with NN-LMS. Just as the gamma correction in image processing, an exponential value 0 < γ < 1 reduces the dynamic range of parameter vector α(n). Each component α i (n) will be closer to 1 no matter it is larger or smaller than 1.

III. ALGORITHM BEHAVIOR MODELING

Let us denote by α * the solution of the unconstrained leastmean-square problem min α E{[y(n) -α x(n)] 2 }. Defining the weight-error vector by v(n) = α(n) -α * , the update equation ( 4) can be written as

v(n + 1) = v(n) + η e(n) D x (n) v(n) + α * γ . (7)
It is assumed that z(n) is stationary, zero-mean and statistically independent of any other signal, which implies that E{z(n) D x (n} = 0. After some calculation and approximation, the expected value of the weight vector error can be expressed recursively by

E{v(n + 1)} = E{v(n)} -η D r (n) R x E{v(n)} -η D s (n) diag{R x E{v(n)} E{v (n)}}. ( 8
)
where D r (n) and D s (n) are diagonal matrices with diagonal entries defined by

r ii (n) = (E{v i (n)}+α * i ) γ -γ (E{v i (n)}+α * i ) γ-1 E{v i (n)} s ii (n) = γ (E{v i (n)} + α * i ) γ-1
We now briefly present the model for the stochastic behavior of the second-order moments of the adaptive weights. As the mean-square estimation error is expressed

E{e 2 (n)} = σ 2 z + trace{R x K(n)}, (9)
we thus derive a recursion for

K(n) = E{v(n) v (n)} start-
ing from the weight error update equation [START_REF] Lawson | Solving Least Squares Problems[END_REF]. Premultiplying equation ( 7) by its transpose, taking the expected value, and using the statistical properties of z(n) leads us to a recursion of the form

K(n + 1) = K(n) -η P 1 (n) K(n) -η K(n) P 1 (n) + η 2 σ 2 z P 2 (n) + η 2 σ 2 z (P 3 (n) + P 3 (n)) + η 2 σ 2 z P 4 (n) -η (P 5 (n) + P 5 (n)) + η 2 P 6 (n) + η 2 (P 7 (n) + P 7 (n)) + η 2 P 8 (n) (10)
where P 1 to P 8 are functions of K(n), E{v(n)} and R x . Due to space constraints in the abstract, their expressions will be given in the final paper.

We now present simulation examples to illustrate the accuracy of the first and second-order models derived in this paper. We consider an example where the correlated inputs x(n) are generated by the AR(1) model x(n) = 1 2 x(n-1)+w(n). The noise w(n) is i.i.d. and drawn from a zero-mean Gaussian distribution so that σ 2

x = 1. Noise z(n) is drawn from a Gaussian distribution with zero mean and σ 2 z = 0.1. The impulse response α * is given by α * = [0.8 0.6 0.5 0.4 0.3 0.2 0.1 -0.1 -0.3 -0.6] Note that it contains negative entries in order to better test the behavior of the algorithm. The mean value E{α i (n)} of each coefficient is shown in Figure 1 and the excess mean square error J emse = trace{R x K(n)} is shown in Figure 2. The simulation curves (blue line) were obtained from Monte Carlo simulation averaged over 100 realizations. The theoretical curves (red line) were obtained from the models. One can notice that all the curves are perfectly superimposed. We then compare NN-LMS and modified NN-LMS with adjusted parameters so that they reach the same steady state error. The results are shown in the Figure 3. Unlike NN-LMS, whose small weights converge more slowly, our modified NN-LMS has more balanced first-order convergence rate. It also has better second-order convergence performance than NN-LMS.

  and α o the solution of the constrained optimization problem. The KKT conditions at the optimum α o can be combined into the following expression[START_REF] Chen | Non-negative least-mean-square algorithm[END_REF],[START_REF] Chen | System identification under non-negativity constraints[END_REF]

Fig. 1 .Fig. 2 .

 12 Fig. 1. Convergence of the coefficients αi(n) for the modified NN-LMS, in the case of correlated input x(n). Two different step sizes are considered: η = 10 -2 on the left, and η = 3 • 10 -3 on the right. The theoretical curves (red line) obtained from the model (8) and simulation curves (blue line) are perfectly superimposed.

Fig. 3 .

 3 Fig. 3. Comparison of NN-LMS and modified NN-LMS. Left: first-order performance. The modified algorithm has a more balanced weight convergence rate. Right: Excess mean square error.


Jie Chen
email: jie.chen@utt.fr

Paul Honeine
email: paul.honeine@utt.fr

A novel kernel-based nonlinear unmixing scheme of hyperspectral images

In hyperspectral images, pixels are mixtures of spectral components associated to pure materials. Although the linear mixture model is the most studied case, nonlinear models have been taken into consideration to overcome some limitations of the linear model. In this paper, nonlinear hyperspectral unmixing problem is studied through kernel-based learning theory. Endmember components at each band are mapped implicitly in a high feature space, in order to address the nonlinear interaction of photons. Experiment results with both synthetic and real images illustrate the effectiveness of the proposed scheme.

I. INTRODUCTION

By providing images with a large number of spectral bands, hyperspectral imagery has been applied with success to a large number of areas, as it allows to obtain more accurate and detailed information than other types of remotely sensed techniques. It has been applied to a number of areas, including environment, land use and agricultural monitoring.

In hyperspectral imagery, each pixel is a mixture of spectral components associated to a number of pure materials, often present in the scene [START_REF] Keshava | Spectral unmixing[END_REF]. The unmixing problem consists of identifying the pure materials, the so-called endmembers, and estimating their contributions (or abundances). The linear mixture model has been widely studied due to its simple physical interpretation and its analytically tractable solution.

Although the linear mixture model has several advantages, there are many situations, involving multiple light scattering effects, which requires nonlinear mixture analysis. Researchers have proposed several unmixing methods to addressr nonlinear mixtures. In [2], [START_REF] Plaza | Nonlinear neural network mixture models for fractional abundance estimation in AVIRIS hyperspectral images[END_REF], the authors designed neural networks for unmixing nonlinearly mixed pixels. However, these methods require their networks to be trained by training data and the quality of training data may affect the performance notably. In [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF], a nonlinear unmixing algorithm based on Bayesian inference was proposed. This method has a high computational complexity and it is dedicated to the bilinear model. In [START_REF] Raksuntorn | Nonlinear Spectral Mixture Analysis for Hyperspectral Imagery in an Unknown Environment[END_REF], the researchers extended the endmember matrix with cross terms of pure signatures to model the light scattering on different materials. However it remains difficult to identify the number and the elements involved in the cross terms.

In the machine learning community, kernel-based algorithms in Hilbert space have been a topic of considerable interest and have gained notable success to deal with nonlinear problems in various applications. Kernel methods have been studied for classification and detection problems for hyperspectral images [START_REF] Camps-Valls | Kernel-based methods for hyperspectral image classification[END_REF], [START_REF] Kwon | Kernel orthogonal subspace projection for hyperspectral signal classification[END_REF]. Kernel-based nonlinear unmixing approaches have also been investigated [START_REF] Broadwater | A comparison of kernel functions for intimate mixture models[END_REF], [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF], [START_REF] Wu | A kernel spatial complexity-based nonlinear unmixing method of hyperspectral imagery[END_REF], by deriving nonlinear algorithms based on linear ones thanks to the kernel trick. Kernels are however applied to the spectral signature of each pure component, as a whole, and independently of interactions between materials, thus operating as nonlinear distortion functions.

In this paper, we formulate a new paradigm where each kernel is applied at a given wavelength, thus combining the contributions of each endmember via a nonlinear mixing rule. This scheme thus has a clear physical interpretation. The abundances can be determined by solving a constrained kernel regression problem. Computer simulations will show the effectiveness of the proposed approach.

II. NONLINEAR MIXTURE MODELS

Suppose that L is the number of spectral bands. Let r = [r 1 , r 2 , . . . , r L ] ⊤ be an L × 1 observed column pixel which is a mixture of R endmembers corrupted by an additive noise. Assume that M is the L × R target endmember matrix denoted by [m 1 , m 2 , . . . , m R ] with each column m i a material signature. For the sake of convenience in the following we also denote the ℓ th row of M, namely, the signatures of the ℓ th wavelength band, by m ⊤ λ ℓ . Let α = [α 1 , α 2 , ..., α R ] ⊤ be a R × 1 abundance column vector associated with the pixel r.

There are muliple situations where light scattering effects can be observed. For instance, multiple scattering effects can occur on complex vegetated surfaces where it is assumed that incident solar radiation is scattered by the scene through multiple bounces involving several endmembers. Bilinear mixture models account for the presence of multiple photon bounces by introducing additional interaction terms in the linear model. The generalized binlinear model studied in [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF] assumes that the observed pixel can be expressed as

r = Mα + R-1 i=1 R j=i+1 γ ij α i α j m i ⊗ m j + n (1)
with n the noise vector and ⊗ the Hadamard product. The following constraints are considered for the parameters

α i ≥ 0, ∀i ∈ i, . . . , R and R i=1 α i = 1 and 0 ≤ γ ij ≤ 1.
If the components of interest are in an intimate association, light interacts with more than one component as it is multiply scattered. Hapke developed a bidirectional reflectance model based on the fundamentals principles of radiative transfer theory [START_REF] Hapke | Bidirectional reflectance spectroscopy, 1, Theory[END_REF]. This model is usually referred to as the intimate mixture model. The reflectance at a given wavelength is approximately given by

r ≈ w 4(µ 0 + µ s ) H(w, µ 0 )H(w, µ s ) ( 2
)
where w is the average single-scattering ambedo, µ 0 is the cosine of the angle of emergence, µ s is the cosine angle of incidence, and H is Chandrasekhar's function for isotropic scattering. The average single-scattering albedo of a homogeneous surface composed by large spherical particles is

w = R i=1 w i G i /(ρ i D i ) R j=1 G j /(ρ j D j )
where G i is the mass fraction for material i, D i is the diameter of the particles of material i, ρ i is the solid density of the particles of material i, and w i is the single-scattering albedo of material i. Their abundances can be estimated from the relative geometric cross section α i = Mi/(ρiDi) j Mj /(ρj Dj ) . Note that the measures must be non-negative and sum-to-one similarly to the constraints in linear mixing model.

III. KERNEL BASED NONLINEAR ANALYSIS OF HYPERSPECTRAL IMAGES

A. Presentation of the method

Considering the linear mixing model where an observed pixel is a linear combination of the endmembers weighted by fractional abundances, each pixel is given by

r = M α + n (3)
We estimate the abundance vector by minimizing the quadratic cost function

J(α) = r -M α 2 = L ℓ=1 r ℓ -α ⊤ m λ ℓ 2 (4)
under the non-negativity and sum-to-one constraints. Now consider a general mixing mechanism specified by

r = F(α, M) + n ( 5
)
where F is an unknown nonlinear function that defines the interaction parameterized by M and α. This requires us to consider a more general problem of the form

ψ * = arg min ψ∈H ψ L ℓ=1 (r ℓ -ψ(α, m λ ℓ )) 2 (6)
where ψ is element of the functional space H ψ . In order to extract the fraction of the pure endmembers, we specify the function ψ as follows

ψ(α, m λ ℓ ) = α ⊤ m λ ℓ + ϕ(m λ ℓ ) ( 7
)
where ϕ is an element of the functional space H ϕ . This means that the mixture mechanism consists of a linear mixing and some kind of nonlinear fluctuation defined by ϕ(•). Let H ϕ be a reproducing kernel Hilbert space (RKHS) of real-valued function on a compact M and •, • be the inner product in this space. We denote κ ϕ : M × M → IR as the kernel function of this space, which satisfies the reproducing property ϕ(m λ ) = ϕ(•), κ ϕ (•, m λ ) for every function ϕ ∈ H ϕ and m λ ∈ M. It can be found that the functional space H ψ of ( 7) is also a RKHS whose kernel function is

κ ψ (m λ ℓ , m λp ) = m ⊤ λ ℓ m λp + κ ϕ (m λ ℓ , m λp ) (8)
By virtue of the Riesz representation theorem, the function ψ(•) of H ψ minimizing (6) can be written as a kernel expansion in terms of available data, namely, the endmember components at each wavelength

ψ(•) = L ℓ=1 β ℓ κ ψ (•, m λ ℓ ) (9)
Replacing this expansion in [START_REF] Camps-Valls | Kernel-based methods for hyperspectral image classification[END_REF], it can be shown the problem becomes linear with respect to the parameters to estimate

β * = arg min β∈IR L r -K ψ β 2 (10)
where K ψ is the Gram matrix whose (ℓ, p)-th element is κ ψ (m λ ℓ , m λp ). With [START_REF] Broadwater | A comparison of kernel functions for intimate mixture models[END_REF], it is can be found that

K ψ = MM ⊤ + K ϕ
where K ϕ is the Gram matrix corresponding to the nonlinear map ϕ. The pixel vector can be thus reconstructed by

r = K ψ β * = M(M ⊤ β * ) + K ϕ β *
Comparing this expression with [START_REF] Kwon | Kernel orthogonal subspace projection for hyperspectral signal classification[END_REF], it can be observed that the second term K ϕ β * approximates the nonlinear effect defined by ϕ(•). The first term, namely, the linear component, allows us to determine the abundance vector α * with

α * = M ⊤ β * (11)
As abundances are determined in this form, taking the nonnegativity and sum-to-one constraints into consideration in this context yields

M ⊤ β ≥ 0 1 ⊤ (M ⊤ β) = 1
where 1 is an R × 1 vector of ones.

The whole optimization problem is now expressed as

min β r -K ψ β 2 (12)
subject to M ⊤ β ≥ 0 (13)

1 ⊤ (M ⊤ β) = 1 (14)
Once β * is calculated, the abundances can be estimated by [START_REF] Hapke | Bidirectional reflectance spectroscopy, 1, Theory[END_REF]. The problem defined by ( 12) to ( 14) is a quadratic programming problem with both equality and inequality constraints. Numerous candidate methods can be found to solve this problem and determine β * , as presented in [START_REF] Chen | Non-negative least-mean-square algorithm[END_REF], [START_REF] Luenberger | Linear and nonlinear programming[END_REF]. In this paper, we do not focus our interests on detailed solution approaches of this optimization problem.

B. Kernel selection and interpretation

The kernel κ ϕ is an important element that should be defined in this approach and it modelizes the nonlinearity in the mixture model. There is no limit in the choice of this kernel function provided that κ ϕ is an inner product in the given Hilbert space. Note that it is not necessary to know explicitly the nonlinear map ϕ. However it is possible to choose a kernel which is closely associated to the application or to the inherent mixing model.

A widely used kernel in machine learning community is the Gaussian kernel κ(m λ ℓ , m λp ) = exp(-m λ ℓm λp 2 /σ 2 ) with bandwidth σ. Some well-known kernels are homogeneous polynomial kernels κ(m λ ℓ , m λp ) = m λ ℓ , m λp q . Polynomial kernels map the input data into a finite dimensional polynomial space. It can be noticed that the second-order polynomial kernel (q = 2) is tightly related to the bilinear mixture model. Moreover, a higher order polynomial implies that multiple scatters between materials is considered.

In what follows, we shall illustrate the relation between selection of κ ϕ and physical models using a simple example. Consider the case where there are three endmembers in the scene and the generalized bilinear mixture model is the mixing mechanism. A pixel is thus expressed by [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF]

r = M α + γ 12 α 1 α 2 m 1 ⊗ m 2 + γ 13 α 1 α 3 m 1 ⊗ m 3 + γ 23 α 2 α 3 m 2 ⊗ m 3 ( 15
)
with constraints α 1 , α 2 and α 3 ≥ 0

3 i=1 α i = 1
where γ 12 , γ 13 and γ 23 are additional attenuation parameters, and ⊗ denotes the Hadamard product. In this case suppose that the second-order homogeneous polynomial kernel is chosen for defining the functional space H ϕ , namely, κ ϕ (m λ ℓ , m λp ) = m λ ℓ , m λp 2 , which can be written in an inner product form

κ ϕ (m λ ℓ , m λp ) = (m 2 λ ℓ,1 , m 2 λ ℓ,2 , m 2 λ ℓ,3 , √ 2m λ ℓ,1 m λ ℓ,2 , √ 2m λ ℓ,1 m λi,3 , √ 2m λ ℓ,2 m λ ℓ,3 ), (m 2 λp,1 , m 2 λp,2 , m 2 λp,3 , √ 2m λp,1 m λp,2 , √ 2m λp,1 m λp,3 , √ 2m λp,2 m λp,3 )
where m λ ℓ,i is the i th component of m λ ℓ . Thus the feature map from input space to the high dimensional space is then described by

Φ : m λ ℓ = (m λ ℓ,1 , m λ ℓ,2 , m λ ℓ,3 ) Φ(m λ ℓ ) = (m 2 λ ℓ,1 , m 2 λ ℓ,2 , m 2 λ ℓ,3 , √ 2m λ ℓ,1 m λv,2 , √ 2m λ ℓ,1 m λ ℓ,3 , √ 2m λ ℓ,2 m λ ℓ,3 ) ∈ IR 6
The basis of this functional space is thus closely related to the bilinear model defined by [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF] except for the additional auto-terms and the coefficient √ 2. From this example it is clear that the kernel κ ϕ (•, •) represents the nonlinear interactions between different material signatures. However, thanks to the kernel, we do not need to map explicitly the original input data into the high dimensional space. This allows us to consider complex photon interaction mechanism just by modifying the kernel κ ϕ . For example, a polynomial kernel with order q = 4 means that more photon reflectance is modeled. Compared with the method in [START_REF] Raksuntorn | Nonlinear Spectral Mixture Analysis for Hyperspectral Imagery in an Unknown Environment[END_REF], which inserts products of pure material signatures as new endmembers, we do not need to extend the endmember matrix by adding such terms. The endmember matrix does not suffer from the problem of choosing cross terms and of large expansion in dimension.

IV. EXPERIMENT RESULTS

A. Experiments with synthetic images

Firstly, we study the proposed algorithm using synthetic data, generated by the bilinear mixture model and the intimate mixture model using three pure materials. We consider three mineral materials with spectral signatures from the USGC library. These spectra consist of 480 contiguous bands ranging from 0.35 to 2.5 µm (See Fig. 1). Two hyperspectral images composed of 2500 pixels were generated respectively with these two models. All the fractions were subject to nonnegativity and constant sum-to-one constraints. Both images were corrupted by additive noise with SNR equal to 20 dB and 30 dB. The Gaussian kernel function was used for the proposed algorithm and preliminary experiments were conducted to finally set the kernel bandwidth σ = 2. For both models, we compared the results of the proposed approach with the following three linear/nonlinear algorithms: 1) Fully constrained least square method in [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF] (denoted by FCLS). This method was derived based on linear mixture model and provides the optimal solution in the linear sense. 2) Nonlinear method in [START_REF] Raksuntorn | Nonlinear Spectral Mixture Analysis for Hyperspectral Imagery in an Unknown Environment[END_REF] (denoted by NLCT). This method extends the endmember matrix with cross terms of pure signatures, to model the light scatter effects.

In the experiment, we added all the second-order cross terms so that it corresponds to the bilinear model. 3) Kernel fully constrained least square method in [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF] (denoted by KFCLS). This is a kernel method which replaces the inner products in the linear mixture problem by kernel functions.

The root mean square error

RMSE = 1 LR jk α ij -αij 2
was used to compare these algorithms. The simulation results of these algorithm are reported in Table I. As the method FCLS was initially derived for linear mixture models, it has a relatively large RMSE compared with the other nonlinear methods. With added second order cross terms, NLCT reduces notably the RMSE for the bilinear mixture model. However, with the intimate model, it marginally improves the performance. This means it does not match this model. KFCLS, which maps signatures of different material into a high dimensional space, performs even worth than the linear method. Compared with all these linear/nonlinear methods, the proposed approach clearly shows a notable improvement of the RMSE.

B. Experiments with a real image

Let us now study these algorithms with a real image, the scene of the Cuprite mining district (Nevada, USA), captured by the airbone visible infrared imaging spectrometer (AVIRIS). A sub-image of size 81 × 101 pixels was chosen here to evaluate the proposed algorithm. The averaged spectral angle between r and reconstructed r was used

SA = 1 N N n=1 θ(r n , rn )
where N is the total number of pixels and θ(r, r) = cos -1 r, r r r

In this experiment, we used VCA algorithm [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF] to extract the endmembers, was R = 3 and R = 5 for illustration purpose.

The results are shown in the Table II.

TABLE II SPECTRAL ANGLE COMPARISONS

R = 3 R = 5 FCLS [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF] 0.0594 0.0233 NLCT [START_REF] Raksuntorn | Nonlinear Spectral Mixture Analysis for Hyperspectral Imagery in an Unknown Environment[END_REF] 0.0591 0.0213 KFCLS [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF] --Proposed 0.0281 0.0183

Note that the results for KFCLS can not be listed as there is no reconstruction model for this kernel method scheme. Compared with the others, the proposed algorithm has the lowest reconstruction errors for 3 and 5 endmembers. The estimated abundances of these four algorithms are illustrated in Fig. 2.

Fig. 1 .

 1 Fig. 1. Three endmembers used in the simulation


Paul Honeine

Cédric Richard

THE ANGULAR KERNEL IN MACHINE LEARNING FOR HYPERSPECTRAL DATA CLASSIFICATION

Keywords: Hyperspectral data, spectral angle, SVM, reproducing kernel, machine learning

Support vector machines have been investigated with success for hyperspectral data classification. In this paper, we propose a new kernel to measure spectral similarity, called the angular kernel. We provide some of its properties, such as its invariance to illumination energy, as well as connection to previous work. Furthermore, we show that the performance of a classifier associated to the angular kernel is comparable to the Gaussian kernel, in the sense of universality. We derive a class of kernels based on the angular kernel, and study the performance on an urban classification task.

INTRODUCTION

Hyperspectral images are now widely available, owing to the development of remote sensing sensors with an improvement in both spectral and spatial resolutions. For instance, airborne sensors provide hyperspectral images with more than a hundred spectral bands and a spatial resolution up to one meter per pixel. Such resolution allows classification of urban structures by virtue of, on the one hand the spatial visualperception, and on the other the spectral physical-features. Finer resolution provides an increase in the dimensionality of the processed data, allowing for a better discrimination between different classes of data, e.g. between trees, roads, bricks, etc. Furthermore, only a limited set of observations with labels is available. However, constructing a classification rule based on a small training set in a high dimensional space is an ill-posed problem.

Kernel-based methods provide the opportunity to overcome these problems, with the Support Vector Machines (SVM) which take advantage of the combination of the regularized structure of the decision rule and the elegant use of a reproducing kernel to measure the similarity between data, independent of their dimension. Recently, SVM were investigated for hyperspectral data classification, and have proven to provide high performance of detection and discrimination. Initially, conventional kernels were used such as the Gaussian kernel [START_REF] Melgani | Classification of hyperspectral remote sensing images with support vector machines[END_REF][START_REF] Rajpoot | SVM optimization for hyperspectral colon tissue cell classification[END_REF], or adjusted to select optimal spectral bands [START_REF] Guo | Customizing kernel functions for svm-based hyperspectral image classification[END_REF], or even combined using spatial and spectral information [START_REF] Camps-Valls | Composite kernels for hyperspectral image classification[END_REF]. Taking into account the spectral signature concept with an invariance to overall energy (e.g. illuminations), the spectral angle [START_REF] Kruse | The spectral image processing system (sips): Interactive visualization and analysis of imaging spectrometer data[END_REF] as a measure of distance has been adapted to operate on a Gaussian kernel in [START_REF] Mercier | Support vector machines for hyperspectral image classification with spectral-based kernels[END_REF][START_REF] Fauvel | Evaluation of kernels for multiclass classification of hyperspectral remote sensing data[END_REF][START_REF] Sap | Spectral angle based kernels for the classification of hyperspectral images using support vector machines[END_REF].

In this paper, we propose the angular kernel, a new measure of similarity between two spectra which is insensitive to their energies. We provide some of its properties and connections to previous work. Performance associated to this kernel are studied in the light of the universality property, comparing it to the consistency well-known kernels such as the classical Gaussian kernel. We derive a class of kernels based on the angular kernel, and illustrate their performance on a real hyperspectral image for classification of urban data. But before, we review the concept of kernel-based machines for hyperspectral data classification.

Kernel-based machines for hyperspectral data

Pioneered by Vapnik's SVM [START_REF] Vapnik | The nature of statistical learning theory[END_REF], kernel-based machines have proven to be successful in many pattern recognition problems. The key issue behind the high generalization ability of SVM is maintained by a complexity control of the solution, tunable by a regularization parameter (often denoted C) which controls the tradeoff between the model simplicity and the fitness to the training data. Furthermore, a central characteristic of these machines is that they can be expressed in terms of inner products of input data. Replacing these inner products with a reproducing kernel provides an efficient way to implicitly map the data into a high, even infinite, dimensional space and apply the original algorithm in this space. Hence, performance depends crucially on the chosen reproducing kernel. By Mercer's theorem, reproducing kernels are positive semi-definite functions, hence can be expressed as an inner product in a high-dimensional feature space. An easy way to construct valid (reproducing) kernels is to apply rules for

i , x j ) = k β k κ k (x i , x j ) R2. Positive Offset κ(x i , x j ) = κ 1 (x i , x j ) + c R3. Product κ(x i , x j ) = k κ k (x i , x j ) R4. Exponential κ(x i , x j ) = exp 1 σ 2 κ k (x i , x j ) R5. Normalization κ(x i , x j ) = κ k (xi,xj ) √ κ k (xi,xi) κ k (xj ,xj )
engineering more complicated kernels from simple ones. Some basis rules are enumerated in Table 1. The first three rules can be combined into the rule: a positive-coefficient polynomial of a reproducing kernel is a valid one. These rules allow to generate most well-known kernels from the linear kernel x i , x j . For instance applying rule R3, or combining rules R2 and R3, with κ k (x i , x j ) = x i , x j gives respectively the homogeneous and inhomogeneous polynomial kernels, while rule R5 provides the normalized linear kernel x i , x j / x i x j . The exponential kernel exp(1 σ 2 x i , x j ) results from rule R4, while the Gaussian kernel exp(-

1 2σ 2 x i -x j 2
) is obtained from normalizing the exponential kernel with rule R5.

In order to take into account the nature of spectral characteristics in the hyperspectral data, the spectral angle [START_REF] Keshava | Distance metrics and band selection in hyperspectral processing with applications to material identification and spectral libraries[END_REF] as a measure of distance is extensively used in the literature, thanks to its invariance to the spectral energy, e.g. illumination. It is defined between two spectra x i and x j as

θ(x i , x j ) = arccos x i , x j x i x j , ( 1
)
where • is the Euclidean distance and •, • its inner product. In order to provide a kernel based on this measure, most work consider it 1 as a distance, and adapt any distance-based kernel for this purpose [START_REF] Mercier | Support vector machines for hyperspectral image classification with spectral-based kernels[END_REF]. The most investigated kernel is the Gaussian kernel [START_REF] Sap | Spectral angle based kernels for the classification of hyperspectral images using support vector machines[END_REF], of the form

κ(x i , x j ) = exp -1 2σ 2 θ(x i , x j ) , (2)
or substituting the angle with its square value in [START_REF] Fauvel | Evaluation of kernels for multiclass classification of hyperspectral remote sensing data[END_REF].

THE ANGULAR KERNEL

Since each spectrum is positive by nature, as well as the ratio in (1), all spectra lie in the positive orthant 2 . We define in this orthant, denoted hereafter by X , the angular kernel as a similarity measure between two spectra, with

α(x i , x j ) = arccos - x i , x j x i x j , (3)
for any pair (x i , x j ) ∈ X 2 . This kernel corresponds to a monotonic increasing transformation of the normalized linear kernel, with values ranging between π/2 and π.

Proposition 1. The angular kernel defined in (3) is a valid reproducing kernel.

Proof. To prove this, recall from trigonometric identities the expansion of the arccos function into an infinite series:

arccos z = π 2 -arcsin z = π 2 - ∞ k=0 (2k)! 2 2k (k!) 2 (2k + 1) z 2k+1 ,
for any |z| ≤ 1. By substituting z with the normalized linear kernel, we obtain the expansion of the angular kernel:

α(x i , x j ) = π 2 + ∞ k=0 (2k)! 2 2k (k!) 2 (2k + 1) x i , x j x i x j 2k+1 (4)
This expansion is a positive-coefficient polynomial of the normalized linear kernel, resulting from rules R1, R2, and R3. Thus, the angular kernel is a valid reproducing kernel.

Therefore, one can use the angular kernel with any kernelbased learning machine, in order to adapt them for (hyper-) spectral data. Next, we give some insights on the geometric structure of the feature space associated to the angular kernel.

Properties of the angular kernel and its feature space

Before proceeding, we establish the connection between the angular kernel and the spectral angle defined in (1), the latter being a distance. For this purpose, recall the trigonometric identity arccos(-u) = πarccos(u). Thus, we have α(x i , x j ) = πθ(x i , x j ).

(

) 5
This equivalence will be useful throughout this paper.

The space associated to the angular kernel has a very rich structure. From the expansion (4), it is obvious that the dimension of the feature space is infinite. Let φ(•) denotes the map induced by this reproducing kernel, mapping the input space to the feature space, i.e. φ : X → H. The norm of the image of any mapped data is

φ(x i ) 2 H = α(x i , x i ) = arccos (-1) = π.
Therefore, all data of the input space are mapped onto the sphere of radius √ π. Moreover, since α(x i , x j ) ≥ 0 on X 2 , the images lie in the positive orthant. The distance is defined as the norm of the difference, with its square value

φ(x i ) -φ(x j ) 2 H = 2 (π -α(x i , x i )) = 2 θ(x i , x i )
, where the last equality follows from [START_REF] Kruse | The spectral image processing system (sips): Interactive visualization and analysis of imaging spectrometer data[END_REF]. In other words, the square distance is equal to twice the spectral angle. Since α(x i , x i ) ∈ [π/2 π], such a distance is upper bounded by √ π. We refer the interested reader to [START_REF] Amari | Improving support vector machine classifiers by modifying kernal functions[END_REF] for a deeper understanding of the geometry of the feature space.

Performance associated to the angular kernel

All these properties provided so far do not give any information about the performance associated to the use of the angular kernel in machine learning, e.g. SVM for classification. Kernels based on the exponential function are the most used ones, such as the Gaussian and exponential kernels. The performance associated to these kernels is often assigned to their expansion in terms of an infinite series of monomials, with fast falling weightings. The angular kernel shares with these kernels such a property, as given in the expansion (4), and thus is likely to give comparable performance.

The generalization abilities of machine learning classifiers of SVM type, independent of the learning scheme, are studied in [START_REF] Steinwart | On the influence of the kernel on the consistency of support vector machines[END_REF], using the concept of universal kernels. The authors show that there exists a certain class of kernels that are consistent for a large class of classification problems, provided a suitably chosen regularization. This class of so-called universal kernels, includes the Gaussian and the exponential kernels. This is formalized here for the angular kernel.

Proposition 2. The angular kernel is a universal kernel on every compact subset of X . Sketch of proof. From ( 4), the angular kernel takes the form

α(x i , x j ) = ∞ k=0 a k x i , x j x i x j k ,
with a k > 0 for all k ≥ 0. Due to [12, Corollary 10], we get the universality of the kernel.

The universality of the angular kernel means that the functions of its associated feature space are capable of approximating all continuous functions on compact subsets in X .

A class of angular kernels

In previous section, we compared the proposed kernel to classical Gaussian and exponential kernels, despite the fact that the angular kernel has no tunable parameter. This property may be advantageous, since we do not need a tuning step to adapt the kernel to the problem under consideration. Next, we provide a class of kernels (with tunable parameters) based on the angular kernel, following the same scheme provided in Section 1.1. The (homogeneous) polynomial kernel associated to the angular kernel takes the form

α p (x i , x j ) = (α(x i , x j )) p ,
for p ∈ IN + . An inhomogeneous counterpart of this kernel can be given by (α(x i , x j ) + c) p , for any positive c. The exponential kernel is defined as the exponential of the angular kernel, up to a multiplicative bandwidth parameter, namely

α e (x i , x j ) = exp 1 σ 2 α(x i , x j ) .
It is worth noting that it is more convenient to use angular values for the bandwidth parameter σ 2 , for instance σ 2 = π which results into a kernel with values within [ √ e e]. In order to construct the equivalent of the Gaussian kernel for the class of angular kernels, we apply the normalization rule R5 to the above exponential kernel. We obtain

α G (x i , x j ) = α e (x i , x j ) α e (x i , x i ) α e (x j , x j ) = exp -1 σ 2 (π -α(x i , x j )) = exp -1 σ 2 θ(x i , x j ))
, where the identity ( 5) is applied. The resulting kernel is equivalent to the Gaussian kernel based on the spectral angle (2), and extensively used in the literature on hyperspectral data. In other words, the angular kernel (3) can be considered as a linear counterpart of the Gaussian kernel (2).

EXPERIMENTAL RESULTS

Data sets are taken with the ROSIS-03 (Reflective Optical System Imaging Spectrometer) provided by the HySenS project. The original hyperspectral image is of the University of Pavia, Italy, with 610-by-340 pixels and 103 frequency bands. For experiments, we took a sub-image of 250-by-250 pixels representing the south-east of the original image, illustrated in Figure 1 (left). Ground truth information about 6 classes are included to train and test the classifiers, as given in Table 2 and illustrated in Figure 1 (middle and right).

For experimentations, an off-the-shelf SVM classifier is used, and applied here in a one-against-all scheme: binary classifiers are trained on each class against the others, while each test observation is assigned to the class with the maximum output. Preliminary experiments were conducted in order to adjust the regularization term in SVM, C, by a search over a logarithmic grid [10 -3 10 3 ] with increment 10 -k . For the angular Gaussian kernel, we need also to adjust the bandwidth parameter, σ 2 , which is determined over a grid search over [π/2 6 π/2], with increment of the form π/2 k .

Experiments were conducted on different kernels from the class of angular kernels. Table 2 summarizes the results associated to three kernels, and give the misclassification error 2.

rates. The polynomial kernel with degrees ranging from 2 to 10 was used, with the cubic one α 3 giving the best performance. Both the exponential and the Gaussian kernels give comparable results, given in the table for the optimal pair of parameters (C, σ 2 ). The angular kernel gives slightly better classification performance for almost all the 6 classes, and an overall better classification rate, even though the other kernels have been tune to their best parameter values.

CONCLUSION

This paper has addressed the problem of classification of hyperspectral data, providing a new class of kernels for machine learning. The analysis has been carried out on the angular kernel, enumerating some of its properties and giving connections to other kernels. Moreover, we showed that this is a universal kernel, resulting into the consistency of the obtained classifier. Preliminary experimental results indicate the adequacy of such reproducing kernels for the classification of hyperspectral data.

Fig. 1 .

 1 Fig. 1. The hyperspectral image (slice at mid spectral-band) considered in this paper (left), with the spatial distribution of the training (middle) and test (right) datasets. The legend is indicated in Table2.

Table 1 .

 1 Some simple rules for engineering a valid kernel from available ones, with β k , c ∈ IR + and σ ∈ IR.

	Rule	Expression
	R1. Linear combination κ(x

Table 2 .

 2 The 6 classes with the ratio of train/test samples, and the misclassification error rates associated to angular kernels.

	Class-name #train #test	α	α 3	α G
	Asphalt	210	2238 7.9 % 6.3 % 11.3 %
	Meadow	188	4259 3.9 % 3.9 %	3.9 %
	Tree	144	889	4.0 % 4.2 %	3.8 %
	Metal sheet 129	647	0.1 % 0.1 %	0.2 %
	Brick	93	741	4.0 % 15 %	8.8 %
	Shadow	12	84	1.0 % 0.5 %	4.7 %
		Overall error: 4.4 % 5.2 % 5.7 %

In fact, they use the absolute value of the spectral angle. However, this quantity is nonnegative for (positive-value) spectral data.

An orthant is the analogue in high dimensional spaces of a quadrant in the plane or an octant in three dimensions.