Sandy Mahfouz

Farah Mourad-Chehade

Paul Honeine

Joumana Farah

Hichem Snoussi

DECENTRALIZED LOCALIZATION USING FINGERPRINTING AND KERNEL METHODS IN WIRELESS SENSOR NETWORKS

Keywords: decentralized processing, fingerprinting, localization, machine learning, wireless sensor networks

Sensors localization has become an essential issue in wireless sensor networks. This paper presents a decentralized localization algorithm that makes use of radio-location fingerprinting and kernel methods. The proposed algorithm consists of dividing the network into several zones, each of which having a calculator capable of emitting, receiving and processing data. By using radio-information of its zone, each calculator constructs, by means of kernel methods, a model estimating the nodes positions. Calculators estimates are then combined together, leading to final position estimates. Compared to centralized methods, this technique is more robust, less energyconsuming, with a lower processing complexity.

INTRODUCTION

Recent advances in electronics and wireless communication technologies have led to the emergence of wireless sensor networks (WSNs). A WSN consists of a large number of tiny autonomous sensors, each of which being capable of sensing, processing, and transmitting environmental information [START_REF] Akyildiz | A survey on sensor networks[END_REF]. One fundamental issue in WSNs is the sensors localization, since all collected data are meaningless if they are not coupled with their correct locations. A trivial solution to this problem consists of equipping each sensor with a positioning device such as the Global Positioning System (GPS) [START_REF] Hofmann-Wellenhof | Global Positioning System: Theory and practice[END_REF]. However, since WSNs are expected to scale up to thousands of sensors, using GPS is impractical. A better solution is to perform relative positioning, where two types of sensors are considered: anchors and nodes. While anchors have known fixed locations, nodes have unknown positions and thus they are localized with respect to anchors.

Many anchor-based localization algorithms have been proposed in literature, especially those based on estimating nodes distances to anchors using the Received Signal Strength Indicator (RSSI) [START_REF] Parker | Vehicular node localization using Received-Signal-Strength Indicator[END_REF][START_REF] Honeine | Localization in sensor networks -a matrix regression approach[END_REF][START_REF] Essoloh | Distributed localization in wireless sensor networks as a pre-image problem in a reproducing kernel hilbert space[END_REF][START_REF] Mourad | Controlled mobility sensor networks for target tracking using ant colony optimization[END_REF][START_REF] Mourad-Chehade | Polar intervalbased localization in mobile sensor networks[END_REF]. RSSI methods have become popular due to their low-power consumption and cost competitiveness. However, using exact RSSIs to estimate distances This work is supported by the Champagne-Ardenne Region in France.

is not always reliable, due to additive noise, obstacles, etc. A more robust RSSI-based approach performs fingerprinting [START_REF] Lin | Performance comparison of indoor positioning techniques based on location fingerprinting in wireless networks[END_REF][START_REF] Robles | 3D fingerprint-based localization for wireless sensor networks[END_REF]. This technique consists in setting a certain number of positions in the environment called offline positions, then building a radio map where each reference position is associated to the RSSI measurements received at this position. A node's position is then determined by processing the RSSI values collected by this node using the constructed radio map. The advantage of location fingerprinting technique is that it takes into account the stationary characteristics of the environment such as multipath propagation, wall attenuation, etc.

In terms of computational concept, algorithms in WSNs can be divided into centralized and decentralized algorithms [START_REF] Nakamura | Information fusion for wireless sensor networks: Methods, models, and classifications[END_REF][START_REF] Sheu | A distributed localization scheme for wireless sensor networks with improved grid-scan and vector-based refinement[END_REF]. The centralized algorithms require the transmission of all measurements to a fusion center (e.g., a sink node) for processing, which often results in unnecessary wasteful energy costs and bandwidth consumption and thereby reduces the lifetime and utility of the network. In decentralized algorithms, local information exchanged between neighboring nodes is processed locally and, hence, information processing is no more limited to one fusion center. Compared to the centralized strategy, the decentralized approach is more robust to failures, since several fusion centers are engaged. It is also less energy consuming and thus more adapted to WSNs limitations. In our previous work [START_REF] Mahfouz | Kernel-based localization using fingerprinting in wireless sensor networks[END_REF], we proposed a centralized localization algorithm combining location fingerprinting and kernel methods. The processing of the collected data and all positions estimations were done by one single entity, i.e., a fusion center.

In this paper, we present a novel decentralized localization algorithm based on our previous work. The proposed algorithm consists in dividing the region of interest into distinct zones depending on the characteristics of the environment, and then setting for each zone a calculator performing computations. The localization algorithm is then run at each calculator, using an offline an online phase. In the offline phase, reference positions are generated in each zone, and RSSI measurements are collected at these positions, leading to a local power map representing the radio of each zone. Then, using kernel methods [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], each calculator defines a local model that associates to the RSSI measurements of its local map the geographical positions where they were collected. In the online phase, each node regularly collects RSSI measurements. Calculators use these RSSIs along with their models to compute local estimates of nodes positions. The local estimates are then combined, leading to one final estimate for each node at a given time. Having multiple calculators for performing estimations highly improves the robustness of the method.The rest of the paper is organized as follows. The proposed method is described in Section 2. In Section 3, the kernel methods are presented. Section 4 illustrates the effectiveness of the method. Finally, Section 5 concludes the paper.

DESCRIPTION OF THE PROPOSED METHOD

The objective of the proposed method is to compute the geographical locations of sensor nodes in an environment of D dimensions, with D = 2 for a 2D environment or D = 3 for a 3D environment. To this aim, two types of sensors are considered: anchors and mobile nodes. Anchors have known fixed locations, denoted by a i , i ∈ {1, . . . , N a }, and mobile nodes are moving with unknown positions, denoted by x j (t), j ∈ {1, ..., N x }, and hence need to be regularly localized. The region of interest is then partitioned into Z distinct zones depending on the characteristics of the environment. Each zone has a calculator, that is a smart device (e.g., a small computer) capable of handling data, performing calculations, and exchanging information with the sensors. Calculators could also be anchors having high computation capabilities. These smart devices have fixed known locations and can be placed anywhere in a zone. Without loss of generality, we consider that zones are rectangular, each having one calculator located at their centers, as shown in Fig. 1. The localization is then performed locally at every calculator.

Each local algorithm consists of two offline and online phases. In the offline phase, n (z) offline positions, denoted by p (z) , ∈ 1, . . . , n (z) , are generated in a uniformly distributed manner in the studied zone z, z ∈ {1, ..., Z}. Then, anchors broadcast signals in the network with the same initial power. Meanwhile, a sensor is temporarily placed at each offline position of each zone, it detects the anchors signals and measures their Received Signal Strength Indicators (RSSIs). All anchors signals are assumed to be received at all offline positions regardless of their zones. Let

ρ (z) = (ρ a1,p (z) , . . . , ρ aN a ,p (z)
) be the RSSI vector of signals received from the N a anchors at the position p (z) , with ∈ {1, ..., n (z) }. Therefore, for each zone z, a local power map is obtained, composed of n (z) couples ρ (z) , p (z) .

Once the local map is set, the algorithm aims at defining a function ψ (z) (•) that associates to each RSSI vector ρ (z) , for ∈ {1, ..., n (z) }, the corresponding position p (z) . Kernel methods provide an elegant framework to determine these functions, as it will be shown in Section 3.

In the online phase, each node travels freely in the net- work and collects meanwhile RSSIs from anchors. Let ρ j (t) be the N a -vector of RSSIs collected by the node j, at time t, from the N a anchors, with j ∈ {1, ..., N x }. Then, each node sends its RSSI vector ρ j (t) to all the calculators located within its communication range. Let I j (t) be the set of indices of the calculators detected by the node j at time t. Calculators referred in I j (t) apply their defined functions ψ (z) (•) to the RSSI vector ρ j (t) to compute local estimates of the node's position as follows:

x(z) j (t) = ψ (z) (ρ j (t)), z ∈ I j (t), j ∈ {1, ..., N x }.
Then, local estimates are sent to the nearest calculator, where a global estimate is given by a combination of all local estimates using:

xj (t) = z∈Ij (t) w z x(z) j (t), j ∈ {1, ..., N x }.
The quantities w z are weights computed with respect to the distances separating the node from the calculators. Indeed, the more the node is close to a given calculator, the more the calculator's estimate is reliable and thus the greater the corresponding weight should be. Since distances to calculators are not available, we will use instead the powers of the messages emitted by the nodes to the calculators while sending their RSSI-vectors. Let ξ z,j (t) be the RSSI of the message sent by the node j to the calculator z at time t. Since signal powers decrease with the increase of their traveled distance, the closer the node j is to the calculator z, the greater ξ z,j (t) theoretically is. Then, the weights w z are given by ξz,j(t) v∈I j (t) ξv,j (t) , z ∈ I j (t). By using several computation points, the proposed method is more robust to failures compared to the centralized scheme where only one fusion center is considered. Indeed, a failure of the fusion center in the centralized method is fatal, whereas it is much less binding in our decentralized method, since many fusion centers, that is calculators, exist. Also, the localization in this method is less energy consuming, since nodes only send information to calculators in their sensing range, whereas in the centralized scheme data must be routed to a single fusion center regardless of the distance to be traveled.

ESTIMATION OF THE LOCAL MODELS

In this section, we take advantage of kernel-based machine learning to define, for each zone z, a model ψ (z) (•) that associates to each RSSI vector ρ (z) , for ∈ {1, ..., n (z) }, the corresponding position p (z) . To this end, let

ψ (z) (•) = ψ (z) 1 (•), . . . , ψ (z)
,D ) . Let us consider a reproducing kernel κ : IR Na × IR Na → IR, and let H be the induced space.

The function ψ

(z) d
∈ H is estimated by minimizing the mean quadratic error between the model's outputs ψ

(z) d (ρ (z) )
and the desired outputs p

(z)
,d , for ∈ {1, ..., n (z) } :

min ψ (z) d ∈H n (z) =1 (p (z) ,d -ψ (z) d (ρ (z) ) 2 + η (z) ψ (z) d 2 H , (1)
where the regularization parameter η (z) controls the tradeoff between the training error and the complexity of the solution.

According to the representer theorem [START_REF] Schölkopf | A generalized representer theorem[END_REF], the optimal function can be written as follows :

ψ (z) d (•) = n (z) =1 α (z) ,d κ(ρ (z) , •). (2)
By injecting ( 2) in (1), we get the dual optimization problem:

min α (z) d p (z) d -K (z) α (z) d 2 + η (z) (α (z) d ) K (z) α (z) d ,
where α

(z) d = α (z) 1,d , . . . , α (z) n (z) ,d , p (z)
d is the vector of entries p (z) ,d for ∈ {1, ..., n (z) }, and K (z) is a n (z) -by-n (z) matrix whose (s, s )-th entry is given by κ(ρ

(z) s , ρ (z)
s ), for s, s ∈ {1, ..., n (z) }. The solution of the above problem is given by the following linear system :

(K (z) + η (z) I n (z) ) α (z) d = p (z) d , (3)
where I n (z) is the n (z) -by-n (z) identity matrix. It is easy to see the impact of the regularization parameter on the wellposedness of the problem, since the matrix between parenthesis is always non-singular for an appropriate value of η (z) . Equation [START_REF] Parker | Vehicular node localization using Received-Signal-Strength Indicator[END_REF] shows that the same matrix (K (z) +η (z) I n (z) ) needs to be inverted in order to estimate each coordinate d. Nevertheless, it is reasonable to collect all D estimations (D being the space's dimension) in a single matrix inversion problem, thus reducing the computational complexity:

α (z) = (K (z) + η (z) I n (z) ) -1 P (z) , (4)
where α (z) = α

D . Using equation ( 2) and the definition of the vector of functions ψ (z) (•), we define now a model for each zone that allows us to find the D coordinates at once as follows:

ψ (z) (•) = n (z) =1 α (z) κ(ρ (z) , •).
We have now defined, for each zone z, a model ψ (z) that associates to each RSSI vector ρ (z) , for ∈ {1, ..., n (z) }, the corresponding position p (z) . The Z models are used for localization as explained at the end of Section 2.

In this paper, we consider the Gaussian kernel given by:

κ(ρ s , ρ s ) = exp -ρ s -ρ s 2 2(σ (z) ) 2 ,
where • is the Euclidian norm and σ (z) the bandwidth of the Gaussian kernel associated to zone/calculator z. This quantity, together with the regularization parameter η (z) , control the degree of smoothness, noise tolerance, and generalization of the solution. Since we have different data and a different model for each zone, σ (z) and η (z) vary from zone to zone. Choosing σ (z) and η (z) for each zone is performed using the cross-validation technique. This approach is a statistical method of evaluating and comparing learning algorithms by dividing data into two segments: one for training the model and the other one for validating it. We employ this technique to choose proper values for η (z) and σ (z) for each model z [START_REF] Stone | Cross-validatory choice and assessment of statistical predictions[END_REF]. In this paper, we use the k-fold cross-validation which is the basic form of cross-validation. The data in each zone z is first partitioned into k roughly equal sized folds. Subsequently, k iterations of training and validation are performed such that, within each iteration, k -1 folds are used for learning and the remaining one for validation. In each iteration, we estimate the error on the validation set for different values of η (z) and σ (z) . In the end, we choose η (z) and σ (z) that give a minimum error for all iterations.

SIMULATIONS

In this section, we evaluate the performance of the proposed method using various scenarios. In the first part, the proposed method is tested with Matlab-generated data. In the second section, the results are shown in the case of real data gathered in a 10m × 10m real indoor environment [START_REF] Zanca | RSSI measurements in indoor wireless sensor networks[END_REF]. In the final section, the results obtained with the proposed method are compared to the ones obtained when performing localization using connectivity information or the Weighted K-Nearest Neighbor (WKNN) algorithm.

Evaluation of the proposed method on Matlabsimulated data

We consider a 120m × 120m area partitioned into four equal zones and we generate 16 static anchors uniformly distributed over the area. For the training phase, 144 offline positions are generated in a uniformly distributed manner over the area, leading to n (z) = 36 offline positions per zone, z = 1, ..., 4.

The RSSI values needed to construct the local power maps are generated using the well-known Okumura-Hata model [START_REF] Medeisis | On the use of the universal Okumura-Hata propagation prediction model in rural areas[END_REF] given by:

ρ ai,p (z) = ρ 0 -10 n P log 10 a i -p (z) + ε (z) i, , where ρ ai,p (z)
(in dBm) is the power received from the anchor a i by the node at position p (z) , ρ 0 is the initial power (in dBm) set to 150, n P is the path-loss exponent set to 4, a ip (z) is the Euclidian distance between the position p (z) of the considered node and the anchor position a i , and

ε (z)
i, is the noise affecting the RSS measures. We use crossvalidation to find the proper values of the tunable parameters η (z) and σ (z) for each model z, considering η

(z) = 2 s with s ∈ {-20, -19, • • • , -1} and σ (z) = 2 s with s ∈ {1, 2, • • • , 10}.
We then consider the localization problem of one mobile node. To this end, we generate a trajectory and calculate the RSSI values collected by the moving node using the same RSSI model as for the offline positions. The powers of messages exchanged between the node and the calculators are also computed using this model and the node is assumed to communicate with all calculators. Fig. 2 shows the generated trajectory and the estimated one using noiseless RSSIs. It is obvious that the estimation is almost exact. The estimation error, measured by the root mean squared distance between the exact positions and the estimated ones, is in this case equal to 0.14m. If we consider the centralized version of this method [START_REF] Mahfouz | Kernel-based localization using fingerprinting in wireless sensor networks[END_REF], the error is equal to 0.16m. To evaluate the robustness of our method against noise, we add a zero-mean Gaussian white noise of standard deviation σ b = 1dB to the RSSI values. The estimated trajectory still follows the actual one and the estimation error is equal to 1.30m, while the centralized version yields 1.23m of errors. Both centralized and decentralized methods have almost the same efficiency when the noise is spread in the same way all over the network. However, besides the robustness of the method, the decentralized method outperforms the centralized one in terms of computation complexity. Indeed, while our proposed method requires only O((n (z) ) 3 ) for the matrix inversion, with n (z) = 36 in this case, the centralized method considers all the offline positions at once with O(n 3 ), where n = 144.

Evaluation of the proposed method using real data

In this section, we study the performance of the new method in the case of real collected data. The set of collected measurements used are available from [START_REF] Zanca | RSSI measurements in indoor wireless sensor networks[END_REF]. The authors deployed 48 uniformly distributed EyesIFX nodes in a room measuring approximately 10m × 10m. Furniture and people in the room cause multi-path interferences affecting the collected RSSI values. We chose 5 nodes to be anchors in our simulation and used the available data to generate 57 new offline positions in order to get a total of 100 offline positions. We partitioned the area into 9 zones and generated one calculator for each zone. To evaluate our algorithm in the case of a moving node, we generated a trajectory using the available data. Fig. 3 shows the estimated trajectories obtained using the proposed method and the centralized one. The estimation errors are equal to 0.28m and 0.37m respectively. By dividing the region into several zones, noises due to the space configuration being higher in some zones than in others would only affect some local models, the others leading to good estimates. This diversity is absent in the centralized model where anywhere noises impact the one and only model.

Comparison with previous techniques

In this section, we compare the results obtained when using the new method to the results obtained when using two different approaches: localization by connectivity and localization using the WKNN algorithm. We first consider the localization by connectivity [START_REF] Shang | Localization from mere connectivity[END_REF] where nodes are localized using only neighboring anchors information. Assume that the transmission range is an ideal disk with equal transmission radius for all anchors. The position of an unknown node detecting its neighboring anchors is given by the intersection of all 1 shows the estimation error (in meters) in the case of noiseless Matlabsimulated data for different numbers of anchors. By adding anchors, the learning process in the proposed method leads to a rougher model with less freedom degrees, which yields increasing estimation errors. However, the error remains less than the connectivity-based one. We consider at second the WKNN algorithm [START_REF] Koyuncu | A 2D positioning system using WSNs in indoor environment[END_REF] using fingerprinting coupled to RSSI-Euclidean distances computations. In WKNN, estimated positions are given by weighted combinations of the positions of the K nearest offline positions. Different weights are used in simulations of WKNN, as listed in Table III of [START_REF] Mahfouz | Kernel-based localization using fingerprinting in wireless sensor networks[END_REF], and the optimal value of K is determined using the cross-validation approach. Table 2 shows the estimation errors (in meters) obtained while estimating the same trajectory as in the case of real data (Fig. 3) with different weights of WKNN and with our method. The simulation results show that our method is more accurate than the WKNN based one. Note that localization by connectivity is not tested here because there are only 5 anchors in the network, which is insufficient for the connectivity technique to operate properly, while our approach is expected to yield good results in such case.

CONCLUSION

In this paper, we propose a novel decentralized localization method using radio-location fingerprinting and kernel methods. The proposed algorithm divides the network into several zones where different local position estimates are computed. A global estimate is then obtained by combining the local ones. Being a decentralized scheme, this method outperforms centralized ones in terms of robustness and computation complexity. Simulation results show that this method is more efficient than the WKNN algorithm for different weight models. It also gives better estimates than localization by connectivity for different number of anchors. Future works will handle improvements of this method, especially in terms of zones definitions and anchors signals ranges.

Fig. 1 .

 1 Fig. 1. Network topology-2 represents the calculators, represents the anchors, and • represents the offline positions.

D

  (•) , where ψ (z) d : IR Na → IR estimates the d-th coordinate. Thus, the function ψ (z) d , for an input ρ (z) , yields the d-th coordinate in p (z) = (p (z) ,1 , . . . , p

Fig. 2 .Fig. 3 .

 23 Fig. 2. Localization in the case of noiseless simulated data.

Table 1 .

 1 Comparison of the proposed method to the connectivity-based one.

	No. of anchors N a	16	49	81	100
	Connectivity	8.43 4.07 3.68 3.25
	Proposed method 0.14 0.20 1.04 1.37

Table 2 .

 2 Comparison of the proposed method to the WKNNbased one, for different weight models.

	"A"	"B" "C"	"D"	"E"	Proposed
	K=2 K=4 K=4 K=10 K=10	method
	0.69 0.63 0.58 0.57	0.60	0.28
	the range disks of these detected anchors. Table


Jie Chen
email: jie.chen@unice.fr

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

Jose-Carlos M Bermudez
email: j.bermudez@ieee.org

Paul Honeine
email: paul.honeine@utt.fr

NON-STATIONARY ANALYSIS OF THE CONVERGENCE OF THE NON-NEGATIVE LEAST-MEAN-SQUARE ALGORITHM

Keywords: Non-negativity constraint, adaptive filtering, non-stationary signal, convergence analysis

Non-negativity is a widely used constraint in parameter estimation procedures due to physical characteristics of systems under investigation. In this paper, we consider an LMS-type algorithm for system identification subject to non-negativity constraints, called Non-Negative Least-Mean-Square algorithm, and its normalized variant. An important contribution of this paper is that we study the stochastic behavior of these algorithms in a non-stationary environment, where the unconstrained solution is characterized by a time-variant mean and is affected by random perturbations. Convergence analysis of these algorithms in a stationary environment can be viewed as a particular case of the convergence model derived in this paper. Simulation results are presented to illustrate the performance of the algorithm and the accuracy of the derived models.

INTRODUCTION

Optimizing a cost function given a set of constraints is a common objective in parameter estimation problems. The constraints are usually imposed by system specifications which provide a priori information on the feasible set of solutions. Estimating parameters subject to constraints poses specific problems in online applications. Nonnegativity is one of the most commonly stated constraints. It is often imposed on the parameters to estimate in order to avoid physically absurd and uninterpretable results. For instance, non-negativity constraints were used for image deblurring [START_REF] Benvenuto | Nonnegative least-squares image deblurring: improved gradient projection approaches[END_REF], deconvolution of system impulse response estimation [START_REF] Lin | Bayesian regularization and nonnegative deconvolution for room impulse response estimation[END_REF] and audio processing [START_REF] Cont | Realtime multiple pitch and multiple-instrument recognition for music signals unsing sparse non-negative constraints[END_REF]. Another similar problem is the non-negative matrix factorization (NMF), which is now a popular dimension reduction technique [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF][START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF][START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF]. This problem is closely related to blind deconvolution, and has found direct application in neuroscience [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multiway data analysis and blind source separation[END_REF] and hyperspectral imaging [START_REF] Berry | Algorithms and applications for approximate nonnegative matrix factorization[END_REF]. Separation of non-negative mixture of non-negative sources has also been considered in [START_REF] Plumbley | Algorithms for nonnegative independent component analysis[END_REF][START_REF] Moussaoui | Separation of non-negative mixture of nonnegative sources using a bayesian approach and MCMC sampling[END_REF].

Over the last fifteen years, a variety of methods have been developed to tackle non-negative least-squares (NNLS) problems. Active set techniques use the fact that if the set of variables which activate constraints is known, then the solution of the constrained leastsquares problem can be obtained by solving an unconstrained one that includes only inactive variables. The active set algorithm of Lawson and Hanson [START_REF] Lawson | Solving Least Squares Problems[END_REF] is a batch resolution technique for NNL-S problems. It has become a standard among the most frequently used methods. In [START_REF] Bro | A fast non-negativity-constrained least squares algorithm[END_REF], Bro and De Jong introduced a modification of the latter, called Fast NNLS, which takes advantage of the special characteristics of iterative algorithms involving repeated use of nonnegativity constraints. Projected gradient algorithms [START_REF] Rosen | The gradient projection method for nonlinear programming. part 1: Linear constraints[END_REF][START_REF] Calamai | Projected gradient methods for linearly constrained problems[END_REF][START_REF] Barzilai | Two-point step size gradient methods[END_REF][START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF] form another class, which is based on successive projections onto the feasible region. In [START_REF] Lin | Projected gradient methods for nonnegative matrix factorization[END_REF], Lin used this kind of algorithm for NMF problems. Low memory requirements and simplicity make these algorithms attractive for large scale problems. Nevertheless, they are characterized by slow convergence rate if not combined with appropriate step size selection. The class of multiplicative algorithms is very popular for dealing with NMF problems [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF][START_REF] Lin | On the convergence of multiplicative update algorithms for nonnegative matrix factorization[END_REF]. Particularly efficient updates were derived in this way for a large number of problems involving non-negativity constraints [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF]. However, these algorithms require batch processing, which is not suitable for online system identification problems.

An LMS-type adaptive algorithm called Non-Negative Least Mean Square (NNLMS), which is suitable for online applications, has been proposed in [START_REF] Chen | Nonnegative least-mean square algorithm[END_REF]. It is based on stochastic gradient descent approach combined with a fixed point iteration scheme. The NNLMS algorithm was studied in [START_REF] Chen | Nonnegative least-mean square algorithm[END_REF] for operation in stationary environments, and was shown to converge to a solution satisfying the Karush-Kuhn-Tucker conditions. In this paper we study the stochastic behavior of the NNLMS algorithm and its normalized counterpart for operation in a nonstationary environment, where the unconstrained solution is characterized by a time-variant mean and is affected by random perturbations. Convergence analysis of these algorithms in a stationary environment can be viewed as a particular case of the convergence model derived in this paper.

NON-NEGATIVE LEAST-MEAN-SQUARE ALGORITHM

Consider the estimation problem where the unknown system is characterized by real-valued observations

y(n) = α * x(n) + z(n), (1)
where (•) denotes the vector transpose, α * = [α * 1 , α * 2 , . . . , α * N ] is the vector of model parameters, and x(n) = [x(n), x(n -1), . . . , x(n -N + 1)] is the input data vector. The input signal x(n) and the additive noise z(n) are assumed stationary and zero-mean.

In certain applications, inherent physical characteristics impose a non-negativity constraint on the estimate α. Therefore, the problem of identifying the optimum non-negative model can be formalized as follows

α o = arg min α J(α) subject to α o i ≥ 0, ∀i, (2)
where J(α) is a continuously differentiable and strictly convex cost function in R N , and α o is the solution to the constrained optimization problem. To solve the problem (2), let us consider its Lagrangian function Q(α, λ) given by [START_REF] Boyd | Convex Optimization[END_REF] Q(α, λ) = J(α)λ α, where λ is the vector of non-negative Lagrange multipliers. The Karush-Kuhn-Tucker conditions must necessarily be satisfied at the optimum defined by α o , λ o , namely,

∇αQ(α o , λ o ) = 0 (3)
α o i [λ o ] i = 0, i = 1, . . . , N (4)
where ∇α stands for the gradient operator with respect to α. Using that ∇αQ(α, λ) = ∇αJ(α)λ, these equations can be combined into the following expression

α o i [-∇αJ(α o )]i = 0, i = 1, . . . , N (5)
where the extra minus sign is just used to make a gradient descent of J(α) apparent. Equations of the form ϕ(u) = 0 can be solved with a fixed-point iteration algorithm by considering the problem u = u + ϕ(u) under some conditions on function ϕ. Implementing this strategy with equation ( 5) leads to the component-wise gradient descent algorithm

αi(n + 1) = αi(n) + ηi(n) αi(n)[-∇αJ(α(n))]i (6)
with ηi(n) a positive step size required to get a contraction scheme and to control the convergence rate. Consider the criterion

JMSE(α) = E{[y(n) -α x(n)] 2
} to be minimized with respect to α so that

α o = arg min α E{[y(n) -α x(n)] 2 } subject to αi ≥ 0, ∀i (7)
The gradient of JMSE(α) with respect to α is ∇αJMSE(α) = 2 (Rxαrxy), where Rx is the autocorrelation matrix of x(n), and rxy is the correlation vector between x(n) and y(n). Following a stochastic gradient approach, the second-order moments Rx and rxy are replaced in (6) by the instantaneous estimates x(n) x (n) and y(n) x(n), respectively. Choosing a positive fixed step size η, defining Dα(n) as the diagonal matrix with diagonal entries given by α(n), leads to the stochastic update given by

α(n + 1) = α(n) + η Dα(n) x(n) e(n) (8)
where the estimation error e(n) = y(n)α (n) x(n), and the algorithm is required to be initialized with all αi(0) > 0. We refer to this method as the Non-Negative LMS algorithm (NNLMS).

Compared with the classical LMS algorithm, the direction of the (instantaneous) gradient vector in ( 8) is modified by the premultiplication by η Dα(n), i.e., the ith component is scaled by η αi(n). The weight vector update is then no longer in the direction of the gradient, as a consequence of the constraints. This effect enables the weight corrections to reduce gradually to zero for weights approaching zero from the positive phase. These weights tend to converge without becoming negative. If a weight that approaches zero turns negative due to the stochastic update, the instant multiplicative weight αi(n) turns the update from the gradient descend strategy to a gradient ascend one, and the weigth tends to turn positive again. Hence, differently from LMS, the NNLMS correction term is a nonlinear function of α(n). This causes the two algorithms to have completely different convergence behavrios A direct extension of the NNLMS is the Normalized NNLMS. Conditionally to α(n), the product e(n)Dx(n) in ( 8) is proportional to the power of the input signal. Hence, setting a constant value for the step size η leads to different weight updates for different signal power levels. This is the same sensitivity to signal power verified in the LMS algorithm. A popular way to address this limitation is to normalize the weight update by the input vector squared 2-norm which yields the Normalized NNLMS update equation

α(n + 1) = α(n) + η x (n) x(n) + Dα(n) α(n) e(n) (9)
A small positive regularization parameter is added to the denominator x (n) x(n) to avoid numerical difficulties.

STOCHASTIC BEHAVIOR STUDY FOR NON-STATIONARY ENVIRONMENTS

The nonlinearity of the weight correction term with respect to the estimated weight vector makes the theoretical analysis of the NNLMS algorithm quite different from that of LMS. Models for the behavior of the algorithm in stationary environments have been studied in [START_REF] Chen | Nonnegative least-mean square algorithm[END_REF]. We now study the stochastic behavior of the NNLMS and its normalized variant with fixed step sizes for a non-stationary environment. To this end, instead of considering a constant system weight vector α * in (1), the system is characterized by a time variant weight vector α * (n) given by

α * (n) = α * o (n) + ξ(n) (10)
with α * o (n) a deterministic time-variant mean trajectory, and ξ(n) a zero-mean random variable with covariance matrix Ξ = σ 2 ξ I that is independent of any other signal. The deterministic trajectory may result from inherent properties of the system, such as cosine or circular behavior [START_REF] Zhao | Diffusion adaptation over networks under imperfect information exchange and nonstationary data[END_REF]. This model leads to a tractable analysis and permits inferences about the behavior of the algorithms in time variant environments, by varying the trajectory α * o (n) and the power σ 2 ξ of ξ(n). For the analyses that follow, we shall define the weight error vector with respect to the mean unconstrained solution α * o (n) as

v(n) = α(n) -α * o (n). (11)
The following analysis is performed for x(n) and z(n) zero-mean stationary Gaussian, and for z(n) white and statistically independent of any other signal. We assume in the subsequent mean weight behavior analysis that the input and weight vectors are statistically independent, according to the well-known Independence Assumption [START_REF] Sayed | Adaptive Filters[END_REF]. We start by studying the behavior of the NNLMS in nonstationary environments. The result will be extended to its normalized variant at the end of this section.

Mean weight behavior analysis

Considering that the estimation error can also be expressed by

e(n) = z(n) -(v(n) -ξ(n))
x(n), and using the relation [START_REF] Lawson | Solving Least Squares Problems[END_REF] in [START_REF] Plumbley | Algorithms for nonnegative independent component analysis[END_REF] the weight error update equation of NNLMS can be written as

v(n + 1) = v(n) + ηz(n)Dx(n)v(n) + ηz(n)Dx(n)α * o (n) -ηDx(n)α * o (n)x (n)v(n) -ηDx(n)v(n)v (n)x(n) + ηDx(n)v(n)ξ (n)x(n) + ηDx(n)α * o (n)ξ (n)x(n) -∆(n). (12)
where Dx(n) is the diagonal matrix with x(n) as diagonal entries, and ∆(n) = α * o (n + 1)α * o (n) is a deterministic vector proportional to the derivative of the mean unconstrained optimal solution. Taking the expected value of [START_REF] Bro | A fast non-negativity-constrained least squares algorithm[END_REF] and noting that the expectations of the second, third, sixth and seventh terms on the r.h.s. are equal to zero by virtue of the natures of z(n) and ξ(n) yields

E{v(n + 1)} = E{v(n)} -ηE Dx(n)α * o (n)x (n)v(n) -ηE Dx(n) v(n) v (n) x(n) -∆(n). ( 13
)
where E {•} denotes the expected value. Using the independence assumption, the second expectation in the r.h.s. of ( 13) can be written as

E Dx(n)α * o (n)x (n)v(n) = Dα * o (n)RxE{v(n)} (14)
The third term is given by

E{Dx(n) v(n) v (n) x(n)} = diag{Rx K(n)} (15)
with K(n) the correlation matrix of the weight error

K(n) = E{v(n) v (n)}.
Hence the mean behavior model ( 13) is expressed by

E{v(n + 1)} = I -η Dα * o (n) Rx E{v(n)} -η diag{Rx K(n)} -∆(n) (16)
This recursion in the variable E{v(n)} requires a model for K(n).

A recursive model will be derived for K(n) in the later subsection. We have observed that a sufficiently accurate and more intuitive mean behavior model can be obtained by neglecting the weight error fluctuation terms and by using the following separation approximation K(n) ≈ E{v(n)} E{v (n)}. A discussion about the validity of this approximation can be found in [START_REF] Chen | Nonnegative least-mean square algorithm[END_REF]. We thus obtain the following first-order model

E{v(n + 1)} = I -η Dα * o (n) Rx E{v(n)} (17)
-η diag Rx E{v(n)}E{v (n)} -∆(n)

Second-order moment analysis

The excess means square estimation error (EMSE) is given by

JEMSE(n) = E (α(n) -α * (n)) x(n)x (n)(α(n) -α * (n))
. Using ( 10)-( 11), the properties of ξ(n), we can write JEMSE(n) as

JEMSE(n) = E (v(n) -ξ(n)) x(n)x (n)(v(n) -ξ(n)) = trace{Rx K(n)} + trace{Rx Ξ} (18)
The term trace{Rx Ξ} is the direct contribution of the random nonstationarity of the system to the EMSE. In order to estimate the EMSE, we need to derive a recursive model for K(n). The same statistic assumptions A1-A4 as in [START_REF] Chen | Nonnegative least-mean square algorithm[END_REF] are used in the following derivation.

Post-multiplying ( 12) by its transpose, taking the expectation, it can be observed that

K(n+1) = K(n)+K ξ (n+1)+K ξ (n+1)+K ∆ (n+1) (19)
Matrix K ξ (n + 1) consists of the expectations of products between the pairs of the first to fifth terms in [START_REF] Bro | A fast non-negativity-constrained least squares algorithm[END_REF], and is not affected by the random fluctuation ξ(n). Matrix K ξ (n + 1) consists of crossproducts where ξ(n) is involved. The last term K ∆ (n) conveys the effect of deterministic variation of the mean of system weights.

Decoupled with the random perturbation, matrix K ξ (n+1) evolution has the same form as the recursive equation of the weight error covariance matrix derived in [START_REF] Chen | Nonnegative least-mean square algorithm[END_REF] for the case of a stationary environment, except that the system weight α * has to be replaced here by its time variant counterpart α * o (n). This leads to the following expression for K ξ (n + 1)

K ξ (n + 1) = η 2 (Dα * o (n) Q(n) Dα * o (n) + Q(n) • K(n)) + ησ 2 z (Dα * o (n)RxDα * o (n) + Rx • K(n)) + (P 1(n) + P 1 (n)) with Q(n) = 2RxK(n)Rx +trace{RxK(n)}Rx, and
P 1(n) = -η (Dα * o (n) RxK(n) + K(n) Rx E{Dv(n)}) + η 2 Dα * o (n) Q(n) E{Dv(n)} + η 2 σ 2 z E{Dv(n)} Rx Dα * o (n),
where o denotes the so-called Hadamard entry-wise product. Notice that the cross-products of the sixth and seventh terms of [START_REF] Bro | A fast non-negativity-constrained least squares algorithm[END_REF] with the other terms lead to zero mean values. This leaves only auto-product terms

K ξ (n + 1) = η 2 (P 2(n) + P 3(n) + P 4(n) + P 4 (n)) (20)
with

P 2(n) = E{ξ (n)x(n)ξ (n)x(n)Dx(n)v(n)v (n)Dx(n)} P 3(n) = E{ξ (n)x(n)ξ (n)x(n)Dx(n)α * o (n)α * o (n)Dx(n)} P 4(n) = E{ξ (n)x(n)ξ (n)x(n)Dx(n)v(n)α * o (n)Dx(n)}
These terms convey the effect of the random part of the environment non-stationarity. Computing the (i, j)-th entry of P 2(n) leads to

[P 2]ij (n) = k l E {ξ k (n)ξ l (n)} E {vi(n)vj(n)} • E {x k (n)x l (n)xi(n)xj(n)} (21)
As E {ξ k (n)ξ l (n)} = 0 only for k = l, using the Gaussian moment factorizing theorem for the term with fourth-order statistics of x(n),

yields k E x 2 k (n)xi(n)xj(n) = Rx trace{Rx} + 2 Rx Rx ij .
We can thus write the result in matrix form

P 2(n) = σ 2 ξ K(n) • (Rx trace{Rx} + 2 Rx Rx) (22)
Similarly, we have the expected values for P 3(n) and P 4(n)

P 3(n) = σ 2 ξ α * o (n) α * o (n) • (Rx trace{Rx} + 2 Rx Rx) P 4(n) = σ 2 ξ E {vN(n)} α * o (n) • (Rx trace{Rx} + 2 Rx Rx)
Finally, it can be easily observed that the last term in [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF] writes

K ∆ (n + 1) = ∆(n)∆ (n) -∆(n) E {v(n + 1)} + ∆(n) - E {v(n + 1)} + ∆(n) ∆ (n) (23)
With these closed form expressions for K ξ (n + 1), K ξ (n + 1), and K ∆ (n+1), we can characterize the transient behavior of EMSE via [START_REF] Lin | On the convergence of multiplicative update algorithms for nonnegative matrix factorization[END_REF]. Taking n → ∞ the steady state performance can also be studied.

Models for Normalized NNLMS

Evaluation of the expected values of the first and second-order update equations for the Normalized NNLMS (9) involves terms containing the denominator x (n)x(n) + . A common approximation that works well for reasonably large N is to neglect the correlation between this term and the others, as it tends to vary much slower [START_REF] Samson | Fixed point error analysis of the normalized ladder algorithms[END_REF][START_REF] Almeida | A statistical analysis of the affine projection algorithm for unity step size and autoregressive inputs[END_REF]. Given this slow variation and usually very small, we approximate x (n) x(n) + by Nσ 2 x , which is also reasonable for large values of N . Therefore, all the models that have been derived above can be used for Normalized NNLMS by replacing the step size η by the equivalent step size η = η Nσ 2

x .

SIMULATION EXAMPLES

In this section, we present simulation examples to illustrate the properties of the algorithms and the accuracy of the derived models. In these examples, the system order is N = 31. The unknown stationary system is defined by α * (st.)

o i = 1 -0.05 i i= 0, . . . , 20 -0.01 (i -18) i = 21, . . . , 30 (24)
The last ten negative coefficients are used to activate the nonnegativity constraint. For the non-stationary case, we consider an unknown response defined by α * (nst.)

o i (n) = α * (st.) o i + |α * (st.) o i | 10 cos 2πn T + 2π(i -1) N + ξi(n)
where the period T of the deterministic sinusoidal component is set to 3500, and ξ(n) is a zero-mean Gaussian random vector with correlation matrix σ 2 ξ I and σ 2 ξ = 0.001. The input x(n) is a correlated signal generated by a first-order AR process defined as follows: x(n) = 0.5 x(n -1) + w(n), with w(n) an i.i.d. zero-mean Gaussian variable. The variance σ 2 w was set to 1 -0.5 2 so that σ 2

x = 1. The noise z(n) is zero-mean i.i.d. Gaussian with variance σ 2 z = 10 -2 . The adaptive weights in αi(0) were all initialized at 10/N for all the realizations. Monte Carlo simulation results were obtained by averaging 100 runs.

Firstly, we ran the NNLMS algorithm and the theoretical models for both stationary and non-stationary environments with the step size η = 0.005. In Figs. 1(a), 1(b), 1(d), 1(e), blue curves show simulation results and red curves show the theoretical predictions. It can be verified that the models accurately predict the algorithm behavior. In Figs. 1(a) and 1(b), all the coefficients satisfy the nonnegativity constraint. Three more theoretical curves have been added to the second-order plots in Figs. 1(d) and 1(e) to illustrate the effect of the random parameter σ 2 ξ . These curves illustrate the expected extra EMSE due to tracking of the variations of the random optimal solution. To preserve visibility of the graph, simulation curves are not presented but it should be noticed that they coincide with the theoretical ones.

Secondly, the Normalized NNLMS algorithm was tested to highlight its properties. Only the results obtained in non-stationary environments are given here. The step size η of the normalized algorithm was set to η = ηN σ 2 x = 0.005 × 30 × 1 = 0.45 in order to make its performance comparable with that of the NNLMS algorithm depicted in Figs. 1(b) and 1(e), where σ 2 x = 1. For the Normalized NNLMS, variance σ 2 w was set in order that σ 2 x = 0.5. The blue curves show simulation results and red curves represent the theoretical predictions. It can be verified that the proposed models accurately predict the algorithm behavior. Moreover, comparison between Figs. 1(b), 1(e) and Figs. 1(c), 1(f) shows that normalization allows the algorithm to basically converge in the same manner, independently of the input power. As previously, three theoretical curves have been added to Fig. 1(f) in order to illustrate the effect of random perturbations on Normalized NNLMS.

CONCLUSION AND PERSPECTIVE

In this paper, we made further inspection into the NNLMS algorithm to address the online system identification problem under nonnegativity constraints. The main contribution concerns the derivation of the first-order and the second-order stochastic behavior models for non-stationary environments. These models include as a particular case the algorithm behavior in a stationary environment [START_REF] Chen | Nonnegative least-mean square algorithm[END_REF]. The accuracy of these models will enable us to use them as important design tools to estimate the transient, steady-state and tracking performances of the proposed algorithms.

Fig. 1 .

 1 Fig. 1. Evolution of system coefficients and excess mean square error. First column: NNLMS in a stationary environment (with additional theoretical curves as a function of ξ(n) in the EMSE figure). Second column: NNLMS in non-stationary environments. Third column: Behavior of Normalized NNLMS in non-stationary environments.


Jie Chen
email: jie.chen@utt.fr

Paul Honeine
email: paul.honeine@utt.fr

NONLINEAR UNMIXING OF HYPERSPECTRAL IMAGES BASED ON MULTI-KERNEL LEARNING

Keywords: Hyperspectral image, nonlinear unmixing, multi-kernel learning

Nonlinear unmixing of hyperspectral images has generated considerable interest among researchers, as it may overcome some inherent limitations of the linear mixing model. In this paper, we formulate the problem of estimating abundances of a nonlinear mixture of hyperspectral data based on a new multi-kernel learning paradigm. Experiments are conducted using both synthetic and real images in order to illustrate the effectiveness of the proposed method.

INTRODUCTION

Spectral unmixing is an important issue to analyze remotely sensed hyperspectral data, due to the relatively low spatial resolution of images. This involves the decomposition of each mixed pixel into a combination of pure endmember spectra, and the estimation of the abundance value for each endmember [START_REF] Keshava | Spectral unmixing[END_REF]. Under the assumption that the endmembers have been identified, hyperspectral image unmixing thus reduces to an inverse problem for estimating the abundances. To be physically interpretable, the latter are required to satisfy two constraints: they must be nonnegative and their sum equal to one.

Although the linear mixture model has been widely used due to its simple physical interpretation and tractable estimation process, there are many situations in which it may not be appropriate and could be advantageously replaced by a nonlinear one. Nonlinear unmixing has generated considerable interest among researchers, and different methods have been proposed to account for nonlinear effects. Using supervised approaches, such as neural network training, is a way to bypass difficulties that come along with unknown complex mixing mechanism [START_REF] Guilfoyle | A quantitative and comparative analysis of linear and nonlinear spectral mixture models using radial basis function neural networks[END_REF][START_REF] Plaza | Nonlinear neural network mixture models for fractional abundance estimation in AVIRIS hyperspectral images[END_REF]. Unsupervised approaches that do not require training samples have also been extensively studied in the literature. In [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF], a nonlinear unmixing algorithm for generalized bilinear mixture model has been proposed based on Bayesian inference. In [START_REF] Raksuntorn | Nonlinear spectral mixture analysis for hyperspectral imagery in an unknown environment[END_REF][START_REF] Nascimento | Nonlinear mixture model for hyperspectral unmixing[END_REF], the authors have extended the collection of endmembers by adding artificial cross-terms of pure spectral signatures to model light scattering effects on different materials. In [START_REF] Nascimento | Unmixing hyperspectral intimate mixtures[END_REF], the authors have addressed the nonlinear unmixing problem with an intimate mixture model.

Nonlinear algorithms operating in reproducing kernel Hilbert spaces (RKHS) have also been investigated. In [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF], the authors have proposed a kernelization of the FCLS algorithm [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF]. Unfortunately, even though this approach considers nonlinear distortions of spectral signatures, it does not explicitly include nonlinear interactions of the endmember spectra. To overcome this drawback, in our previous work [START_REF] Chen | A novel kernel-based nonlinear unmixing scheme of hyperspectral images[END_REF], we have proposed a new framework based on reproducing kernel machinery in order to possibly take nonlinear photon interactions into account. The family of mixing models that was considered in that paper are partially linear models, which combines a linear mixing model parameterized by the abundances to be estimated, with a nonlinear nonparametric function in a RKHS. In this paper, we extend this approach by formulating the problem as a multi-kernel learning paradigm, in order to optimize the tradeoff between the linear and nonlinear terms in the mixing model. Experiments are conducted using both synthetic and real images in order to illustrate the effectiveness of the proposed method.

PRESENTATION OF THE ALGORITHM

Notations

Suppose that in a given hyperspectral image, there exist R significant endmembers with spectral signature m i ∈ R L , where L denotes the number of spectral bands. Note that, usually, we have R ≪ L. Let r ∈ R L be an observed hyperspectral pixel, and α ∈ R R the unknown abundance fractions associated to the latter. Let M = [m 1 , . . . , m R ] ∈ R L×R be the matrix of the endmembers. For the sake of convenience, the ℓ-th row of M is denoted by m ⊤ λ ℓ ∈ R L , that is, the vector of the endmember signatures at the ℓ-th wavelength band.

Formulation of the nonlinear unmixing problem

There are many situations where it is inappropriate to assume a linear model for hyperspectral data unmixing. See, e.g., the intimate mixture model [START_REF] Keshava | Spectral unmixing[END_REF]. This motivate us to consider general mixing mechanisms of the form

r = Ψ(α, M ) + n (1)
with Ψ an unknown nonlinear function that defines the interactions between the endmembers in matrix M , and n a measurement noise vector. To address the estimation of the abundance vector α, we suggest to consider the following inverse problem

ψ * = arg min ψ∈H 1 2 ψ 2 H + 1 2µ L ℓ=1 (r ℓ -ψ(m λ ℓ )) 2 (2)
with H a given functional space, and µ a positive parameter that controls the trade-off between regularization and fitting.

Clearly, this approach may fail if the functionals ψ of H cannot be adequately and finitely parameterized. Kernel-based methods can lead to efficient and accurate resolution of the inverse problem [START_REF] Guilfoyle | A quantitative and comparative analysis of linear and nonlinear spectral mixture models using radial basis function neural networks[END_REF]. In order to extract the mixing ratios of the endmembers, we define the function ψ in (2) as a partially linear model, which combines a linear mixing model parameterized by the abundances to be estimated, with a nonlinear nonparametric function, namely,

ψ(m λ ℓ ) = α ⊤ m λ ℓ + ψ nlin (m λ ℓ )
subject to α 0 and

1 ⊤ α = 1 (3)
with ψ nlin any real-valued function on a compact M, found from an RKHS denoted by H nlin . Let κ nlin be its reproducing kernel. The constraints in (3) ensure that the abundance fractions are nonnegative, and that their sum is equal to one. We propose to conduct hyperspectral data unmixing by solving the following optimization problem

ψ * , u * = arg min ψ,u 1 2 1 u ψ lin 2 Hlin + 1 1 -u ψ nlin 2 Hnlin + 1 2µ L ℓ=1 e 2 ℓ where ψ = ψ lin + ψ nlin with ψ lin (m λ ℓ ) = h ⊤ m λ ℓ subject to e ℓ = r ℓ -ψ(m λ ℓ ) and h 0 0 ≤ u ≤ 1 (4)
where u is used to optimize the tradeoff between the linear model ψ lin and the nonlinear nonparametric function ψ nlin . Learning simultaneously the parameter u, and functions ψ lin and ψ nlin , in a single optimization problem is known as the multi-kernel learning problem. See [START_REF] Rakotomamonjy | Sim-pleMKL[END_REF] and references therein. The rest of this section is devoted to the formulation and resolution of a convex optimization problem. Comparing equations ( 3) and ( 4), it is important to note that the sumto-one constraint 1 ⊤ h = 1 has been given up due to the antagonistic action of the parameter u. To get α after reaching convergence, an explicit normalization of h by the sum of its entries can be performed.

Convexity of the problem

Problem ( 4) is a convex optimization problem by virtue of the convexity of the function f (u, ψ) = ψ 2 H /u over H × R * + . See [START_REF] Boyd | Convex Optimization[END_REF]. This allows us to formulate the following two-stage optimization problem, with respect to ψ and u successively, in order to solve problem (4)

min u J(u) subject to 0 ≤ u ≤ 1 (5)
where J(u) is defined by [START_REF] Nascimento | Nonlinear mixture model for hyperspectral unmixing[END_REF], see next page. The function J(u) in subproblem ( 6) is defined as the optimum of a family of convex functions in (u, ψ), subject to convex constraints over ψ. As proven in [START_REF] Boyd | Convex Optimization[END_REF], it turns out that J(u) is convex in u and, as a consequence, that the constrained optimization problem (5) is convex.

Solution to the problem

By the strong duality property, we can derive a dual problem that has the same solution J(u) = F (u, ψ * ) as the primal problem [START_REF] Nascimento | Nonlinear mixture model for hyperspectral unmixing[END_REF]. Let us introduce the Lagrange multipliers β ℓ and γ r . The Lagrange function associated with problem ( 6) can be written as

G = 1 2 1 u h 2 + 1 1 -u ψ nlin 2 Hnlin + 1 2µ L ℓ=1 e 2 ℓ - L ℓ=1 β ℓ (e ℓ -r ℓ + ψ(m λ ℓ )) - R r=1 γ r h r (8)
with γ r ≥ 0, where we have used ψ lin 2 Hlin = h 2 . The conditions for optimality of G with respect to the primal variables are given by

     h * = u L ℓ=1 β * ℓ m λ ℓ + γ * ψ * nlin = (1 -u) L ℓ=1 κ nlin (•, m λ ℓ ) β * ℓ e * ℓ = µ β * ℓ (9)
As expected, note from equation ( 9) that u and 1 -u control proportions of the linear and nonlinear components in the model. By substituting equation ( 9) into (8), we get the dual problem [START_REF] Rakotomamonjy | Sim-pleMKL[END_REF], see next page, where K nlin is the Gram matrix associated with ψ nlin , defined by

[K nlin ] ℓ,p = κ nlin (m λ ℓ , m λp ). (10)
Problem ( 11) is a quadratic program (QP). Numerous candidate methods exist to solve it, such as interior point, active set and projected gradient. Function J(u) can thus be evaluated at any point u 0 by solving this QP problem, with u = u 0 .

Once the problem [START_REF] Rakotomamonjy | Sim-pleMKL[END_REF] has been solved, the variable u in the problem (5) is updated using a gradient descent strategy. The existence and computation of the derivatives of supremum functions such as J(u) have been largely discussed in

J(u) =            min ψ F (u, ψ) = 1 2 1 u ψ lin 2 Hlin + 1 1 -u ψ nlin 2 Hnlin + 1 2µ L ℓ=1 e 2 ℓ
subject to e ℓ = r ℓ -ψ(m λ ℓ ) with ψ = ψ lin + ψ nlin and ψ lin (m λ ℓ ) = h ⊤ m λ ℓ with h 0 (6)

J(u) =            max β,γ G ′ (u, β, γ) = - 1 2 β γ ⊤ K u + µI uM uM ⊤ uI β γ + r 0 ⊤ β γ subject to γ 0 with K u = uM M ⊤ + (1 -u)K nlin ( 11
)
the literature [START_REF] Bonnans | Optimization problems with perturbations: A guided tour[END_REF]. The derivatives of J at u 0 can be calculated as if the corresponding minimizer (β * 0 , γ * 0 ) was independent of u 0 , i.e.,

dJ(u) du u=u0 = ∂G ′ (u, β * 0 , γ * 0 ) ∂u u=u0 = - 1 2 M ⊤ β * 0 + γ * 0 2 -β * 0 ⊤ K nlin β * 0 ( 12
)
An appropriate selection of the step size allows the gradient descent algorithm to ensure that the constraint 0 ≤ u ≤ 1 is satisfied. After a given stopping criterion is met, such as small modulus of the derivative [START_REF] Boyd | Convex Optimization[END_REF], pixel reconstruction can be performed using

r * = [ψ * (m λ1 ), . . . , ψ * (m λL )] ⊤ (13)
with ψ * (m λ ℓ ) = m ⊤ λ ℓ h * + ψ * nlin (m λ ℓ ) defined in [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF]. Finally, the estimated abundance vector is given by

α * = M ⊤ β * + γ * 1 ⊤ (M ⊤ β * + γ * ) . (14)

EXPERIMENTS

We shall now conduct some simulations to validate the proposed unmixing algorithm, and to compare it with state-ofthe-art methods, using both synthetic and real images.

Experiments on synthetic images

Let us first report some experimental results on synthetic images. The materials considered in the following scenes are epidote, kaolinite, buddingtonite, alunite, calcite, almandine, jarosite and lepidolite, whose spectra consist of 420 contiguous bands. In the first scene, only the first three materials were selected to generate images. In the second scene, the first five materials were used. In the third scene, the eight materials were considered. For each scene, three 50-by-50 images were generated with different mixture models, each providing N = 2500 pixels for evaluating and comparing performances. These three models were the linear model, the generalized bilinear mixture model with attenuation factors γ ij equal to 1, and a post-nonlinear mixing model (PNMM) defined by r = (M α) 0.7 + n, where ( • ) 0.7 denotes the exponential value 0.7 applied to each entry of the vector. All these images were corrupted with an additive white Gaussian noise n with SNR of 30 dB, except in the case illustrated by Figure 1 that is described hereafter. For our algorithm, the Gaussian kernel

κ nlin (m λp , m λ ℓ ) = exp(-m λp -m λ ℓ 2 /2σ 2 )
was selected to be the reproducing kernel of H nlin .

For illustration purpose only, Figure 1 depicts the evolution of u and 1 -u that was observed during a run of the twostage optimization algorithm derived from ( 5)-( 6). These preliminary simulations were conducted with a noiseless mixed pixel composed of 5 endmembers. The step sizes were adjusted in order to provide a clear demonstration. It can be observed that, in the case of the linear mixing, the weighting coefficient u converged towards 1. This means that, as expected, the ummixing algorithm only captured the linear component with ψ lin . It converged to other values when confronted with the bilinear and post-nonlinear mixing models.

Our approach was compared with the fully constrained least square method (FCLS) [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF], the kernel fully constrained least square method (KFCLS) [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF], the extended endmembermatrix method (ExtM) [START_REF] Raksuntorn | Nonlinear spectral mixture analysis for hyperspectral imagery in an unknown environment[END_REF], and the Bayesian algorithm derived for generalized bilinear model (BilBay) [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF]. The root mean square error was used to compare these five algorithms. Preliminary runs were performed using independent data to tune their parameters. The results are reported in Table 1. Our method was the most efficient one in these experiments.

Experiments with AVIRIS images

Let us now illustrate the performance of the proposed algorithm, when applied to a well known image captured on the Cuprite mining district (NV, USA) by AVIRIS. A sub-image of 30×101 pixels was chosen to evaluate and compare the performance of the algorithms that have been previously enumerated. This area of interest has L = 189 spectral bands. VCA algorithm was used to extract four endmembers. To evaluate the performance, the averaged spectral angle between the original r and the reconstructed r * pixel vectors was used

SA = 1 N N n=1 θ(r n , r * n )
with N is the number of pixels. The performance of each approach is reported in Table 2. Note that KFCLS was not considered in these tests as there is no possible direct reconstruction of pixels. Clearly, in this example, our algorithm had lower reconstruction errors than the other approaches.

Fig. 1 .

 1 Fig. 1. Evolution of u (curves with 'o') and 1 -u (curves with 'x') for the three models. From left to right: linear, bilinear, PNMM

Table 1 .

 1 RMSE comparison

			R = 3			R = 5			R = 8
		linear	bilinear PNMM	linear	bilinear PNMM	linear	bilinear PNMM
	FCLS	0.0037	0.0758	0.0604 0.0134	0.1137	0.1427 0.0148	0.0930	0.1079
	ExtM	0.0079	0.0312	0.0601 0.0157	0.0575	0.1427 0.0173	0.0560	0.1126
	KFCLS	0.0054	0.2711	0.2371 0.0200	0.2051	0.1955 0.0216	0.1431	0.1274
	BilBay	0.0384	0.0285	0.1158 0.0585	0.0441	0.1741 0.0448	0.0369	0.1159
	Proposed	0.0104	0.0315	0.0230 0.0196	0.0288	0.0346 0.0185	0.0221	0.0291

Table 2 .

 2 Spectral angles comparison

		FCLS	ExtM	BilBay Proposed
	SA	0.0206 0.0126 0.0154	0.0110


Nguyen Hoang Nguyen
email: nguyenh@oca.eu

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

Paul Honeine
email: paul.honeine@utt.fr

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

HYPERSPECTRAL IMAGE UNMIXING USING MANIFOLD LEARNING METHODS DERIVATIONS AND COMPARATIVE TESTS

In hyperspectral image analysis, pixels are mixtures of spectral components associated to pure materials. Although the linear mixture model is the mostly studied case, nonlinear techniques have been proposed to overcome its limitations. In this paper, a manifold learning approach is used as a dimensionality-reduction step to deal with non-linearities beforehand, or is integrated directly in the endmember extraction and abundance estimation steps using geodesic distances. Simulation results show that these methods improve the precision of estimation in severely nonlinear cases.

INTRODUCTION

In hyperspectral imagery, a pixel is a mixture of spectral components associated with a number of pure materials present in the scene. In order to reveal underlying information, one needs to identify the endmembers present in the measured mixed pixels and derive the relative proportions of these materials. Under the assumption that the endmembers have been determined a priori using some off-the-shelf extraction techniques, see [START_REF] Winter | N-FINDR: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF][START_REF] Chang | A new growing method for simplex-based endmember extraction algorithm[END_REF][START_REF] José | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF][START_REF] Boardman | Geometric mixture analysis of imaging spectrometry data[END_REF], unmixing of hyperspectral images then consists of estimating their abundances [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF]. In [START_REF] Honeine | Geometric unmixing of large hyperspectral images: A barycentric coordinate approach[END_REF], these two problems are combined using a geometric formulation. The abundances are defined as barycentric coordinates in a reduced-dimensional space, and can be derived using either a ratio of simplex volumes, or a ratio of distances.

Most existing approaches are focused on linear models due to their inherent simplicity from conceptual and implementational points of view. In real images, the interaction between materials can cause nonlinear effects that may influence the precision in abundance calculation, which can violate the sum-to-one constraint and the non-negativity constraint. A possible strategy is, first, to use nonlinear methods for dimensionality reduction, such as ISOMAP [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] and LLE [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF], and then to apply usual linear methods on reduced data. An interesting example is the combined use of geodesic distances and N-FINDR in [START_REF] Heylen | Non-linear spectral unmixing by geodesic simplex volume maximization[END_REF].

In this paper, we generalize the approach described in [START_REF] Heylen | Non-linear spectral unmixing by geodesic simplex volume maximization[END_REF] thanks to the use of geodesic distances with the geometric formulations for abundance estimation proposed in [START_REF] Honeine | Geometric unmixing of large hyperspectral images: A barycentric coordinate approach[END_REF]. This paper is organized as follows. Section 2 is dedicated to endmember extraction and abundance estimation using a geometric approach. Section 3 introduces a new framework based on manifold learning to perform nonlinear unmixing, i.e., endmember extraction and abundance estimation. Section 4 provides some experimental results. Some concluding remarks are presented in Section 5.

ENDMEMBER EXTRACTION AND ABUNDANCE

ESTIMATION: A GEOMETRIC FRAMEWORK

Basic linear mixing and abundance estimation

The basic linear mixing model can be expressed as

r = M α + n (1)
with n a noise vector, r = [r 1 , . . . , r L ] an L-by-1 column pixel supposed to be a mixture of R endmember spectra. Each column m i of the L-by-R matrix M is the spectrum of an endmember, and α is the R-by-1 vector of abundances.

In practice, the spectral dimension is reduced to (R -1) in order that the following mixing model be satisfied

r = R i=1 α i m i (2)
subject to the equality constraint R i=1 α i = 1 and inequality constraints α i ≥ 0. Suppose that the endmembers m i have been extracted using any off-the-shelf technique. We shall now give a direct geometric framework to estimate the abundances in the same low-dimensional space. See [START_REF] Honeine | Geometric unmixing of large hyperspectral images: A barycentric coordinate approach[END_REF] for extensive details. Estimating α in (2) subject to the sum-to-one constraint can be performed by solving the following system

1 1 . . . 1 m 1 m 2 . . . m R α = 1 r . (3)
Using Cramer's rule, the solution can be written as

α i = det 1 . . . 1 . . . 1 m 1 . . . r . . . m R det 1 . . . 1 . . . 1 m 1 . . . m i . . . m R , (4)
for all i = 1, . . . , R. This results can also be expressed by geometric formulation. Consider the simplex S = {m 1 , m 2 , . . . , m R }. Its oriented volume can be calculated as follows

V S = 1 (R -1)! det 1 1 . . . 1 m 1 m 2 . . . m R (5)
Let \ denotes the set difference operator, S \ {m i } ∪ {r} is the set S with m i replaced by r. The calculation of each abundance coefficient α i can thus be performed as

α i = V S \{mi}∪{r} V S (6)
for all i. This can also be written in term of distances

V S = 1 (R -1)! δ(m i )V S \{mi} (7)
where δ(m i ) is the signed distance between the vertex m i and the subspace spanned by the other vertices of S. This leads us to

α i = δ(r) δ(m i ) (8)
for all i. As shown in [START_REF] Honeine | Geometric unmixing of large hyperspectral images: A barycentric coordinate approach[END_REF], by using these equations, one can easily integrate abundance estimation into endmember extraction algorithms [START_REF] Winter | N-FINDR: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF][START_REF] Chang | A new growing method for simplex-based endmember extraction algorithm[END_REF][START_REF] José | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF][START_REF] Boardman | Geometric mixture analysis of imaging spectrometry data[END_REF].

On non-negativity constraint violation

The non-negativity constraint is violated if the pixel under consideration lies outside the simplex S formed by endmembers. In equations ( 6)-( 8), volumes and distances are oriented. Their sign can be either positive or negative, depending on the order of the sequence defined by simplex vertices in the case of volume-based expressions, or on which side of the hyperplane the pixel lies in the case of distance-based expressions. Thus, equations ( 6)-( 8) are appropriate for pixels that lie inside the simplex S. A negative α i means that the pixel is outside S. An alternative approach, which lies beyond the scope of this paper, should be used to process these pixels.

MANIFOLD LEARNING FOR ENDMEMBER EXTRACTION AND ABUNDANCE ESTIMATION

The hypothesis model used by endmembers extraction methods and abundance estimation processes, such as those described above, are basically linear. Due to nonlinear physical effects such as multiple light scattering, this assumption may be uncorrect. Manifold learning can be used to overcome this drawback, by mapping the data into a lower-dimension manifold and returning the desired simplex. Manifold learning techniques considered hereafter are the basic ISOMAP [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] and LLE [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF]. We shall now present some endmember extraction methods combined with abundance estimation that are based on manifold learning.

N-FINDR

The N-FINDR algorithm aims to find the endmembers by maximizing the volume of the simplex formed by pixels. The apices that define this simplex are the endmembers, provided that there exists a pure pixel for each endmember in the data set. A nonlinear dimensionality reduction is first performed with LLE or ISOMAP. The former maps data into a single global coordinate system of lower dimensionality, and can thus be used as a preprocessing step for endmember extraction and abundance estimation in this new coordinate system. As shown below, the latter can also be jointly used with equations ( 6) and ( 8) for endmember extraction and abundance estimation as it returns geodesic distances. Volume calculation is usually accomplished by computing the determinant of a simplex matrix. It can also be written in terms of the Cayley-Menger determinant as shown independently in [START_REF] Honeine | Geometric unmixing of large hyperspectral images: A barycentric coordinate approach[END_REF][START_REF] Heylen | Non-linear spectral unmixing by geodesic simplex volume maximization[END_REF]. Let d ij be the Euclidean distance between endmember spectra m i and m j . The square volume of the simplex defined by the m i 's is given by

V 2 S = (-1) R 2 R-1 ((R -1)!) 2 det(C 1,...,R ) (9)
with

C 1,...,R = D 2 1,...,R 1 1 0 (10)
where D 2 1,...,R is the squared distance matrix whose (i, j)-th entry is D 2 R (i, j) = d 2 ij for all i, j = 1, . . . , R. Note that we can also rewrite [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] as

det(C 1,...,R ) = -(d 1 C -1 2,...,R d 1 ) det(C 2,...,R ) with d 1 = (d 2 12 , . . . , d 2 1R , 1
). Finally, we can calculate the distance between the vertex m 1 and the subspace spanned by the other vertices (m 2 , . . . , m R ) as follows

δ(m 1 ) = d 1 C -1 2,...,R d 1 2 1 2 (11)
Distances in equation ( 9) can be either Euclidean distances, or geodesic distances calculated with ISOMAP algorithm. This allows us to calculate α via ( 6) or ( 8) based on manifold learning considerations.

SGA

The SGA algorithm successively determines the endmembers as follows. First, the algorithm starts with a random point t, and then searchs for another point that maximizes the volume of the simplex S 1i = {t, r i }. This point is chosen as the first endmember and labeled m 1 . SGA continues with the simplexes S 2i = (m 1 , r i ). The point r i that maximizes the volume of this simplex is selected as the second endmember and labeled m 2 . The searching process continues until the R required endmembers are determined.

Using the Cayley-Menger determinant with geodesic distances calculated by ISOMAP leads to an improved SGA based on manifold learning. This algorithm jointly calculates volumes by using [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF], and the abundances by using [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF]. LLE algorithm should also be used, but as a preprocessing stage for dimensionality reduction since it does not return geodesic distances. In this case, equations ( 6) and ( 8) can used for abundance estimation in this new coordinate system.

VCA

The VCA algorithm is based on the fact that the affine transformation of a simplex is also a simplex. It then projects the data onto a direction that is orthogonal to the subspace spanned by the endmembers that were found previously. The farthest point gives the new endmember. The algorithm iterates until it reaches the desired number of endmembers. Basically, VCA processes the distance between a point and the subspace spanned by the endmembers. Again, we can use geodesic distances returned by ISOMAP with ( 8) and (11) to get a manifold-learning based algorithm that jointly determines the endmembers and estimates the abundances.

SIMULATION RESULTS

We shall now compare all these algorithms using an artificial data set that consists of random samples in a two-dimensional simplex, transformed into a three-dimensional nonlinear manifold by projection on a swissroll structure. See [START_REF] Heylen | Non-linear spectral unmixing by geodesic simplex volume maximization[END_REF] for details on the experimental setup. The non-linearity of the swissroll data is parameterized by a variable σ. The endmembers are also in the dataset. The coordinate of a data point r i as a function of the abundance local α i can thus be expressed by

r i1 = α i1 sin(σα i1 ) + 1 r i2 = α i1 cos(σα i1 ) + 1 r i3 = α i2 + 1 (12)
Following the sum-to-one constraint, the abundance of the third endmembers can be generated by

α i3 = 1 -(α i1 + α i2 ).
By setting a single abundance to one, the others to zero, one can find the endmembers m 1 = (sin(σ) + 1, cos(σ) + 1, 1)

m 2 = (1, 1, 2) m 3 = (1, 1, 1) (13)
Swissroll data unmixing was performed over a 1000-sample set using the conventional and the manifold-learning based algorithms. The number of endmembers was set to R = 3. Endmember extraction was performed with N-FINDR, SGA and VCA, in a linear and a nonlinear way. Linear processing consisted of performing first a PCA to map the data into a two-dimensional space, and then applying one of the above three algorithms. LLE and ISOMAP were considered for nonlinear processing. The former was used as a preprocessing step, prior to N-FINDR, SGA and VCA, which were jointly applied with ( 6) and ( 8) to extract endmembers and evaluate their fractions of abundance in mixed pixels. The latter was applied to compute the geodesic distances in ( 9)-( 11), which were used by N-FINDR, SGA and VCA to extract the endmember, and to evaluate the abundances via ( 6) and ( 8).

These algorithms were compared for σ values in [0, 5]. The upper graphic in Figures 123shows the mean of the minimum spectral angle between estimated endmembers m i and the ground truth m i .

θ = 1 R R i=1 min j arccos m i • m j m i m j (14)
It can be observed that the nonlinear algorithms yield acceptable error values over the entire range of σ, and perform better than the linear algorithms. As the nonlinearity parameterized by σ increases, they however fail to retrieve all the endmembers. The lower graphic in Figures 123shows the error between the original and the estimated abundances defined by

E = 1 R R k=1 min k 1 N N i=1 | α ik -α ik | (15)
where N = 1000 is the size of the data set. Again, nonlinear algorithms have better performance than their linear counterparts. Figures 45present the same results in a different way, which allows to compare the performance of N-FINDR, SGA and VCA algorithms when associated with ISOMAP and LLE, respectively. It appears that ISOMAP performs better than LLE, but needs more time-consuming calculations.

CONCLUDING REMARKS

In this paper, we studied several manifold-learning based approaches for joint endmember extraction and abundance estimation in hyperspectral images. Simulation showed that they perform better than usual linear approaches in nonlinear unmixing problems. These results were confirmed by experiments on well-known hyperspectral scenes, which will be presented in the camera-ready paper.

Fig. 1 :Fig. 2 :Fig. 3 :Fig. 4 :Fig. 5 :

 12345 Fig. 1: N-FINR algorithm with PCA (linear), ISOMAP and LLE


Maya Kallas

Paul Honeine

Clovis Francis

Hassan Amoud

PREDICTION OF TIME SERIES USING YULE-WALKER EQUATIONS WITH KERNELS

Keywords: autoregressive model, Yule-Walker equations, expected kernels, pre-image problem, nonlinear model

The autoregressive (AR) model is a well-known technique to analyze time series. The Yule-Walker equations provide a straightforward connection between the AR model parameters and the covariance function of the process. In this paper, we propose a nonlinear extension of the AR model using kernel machines. To this end, we explore the Yule-Walker equations in the feature space, and show that the model parameters can be estimated using the concept of expected kernels. Finally, in order to predict once the model identified, we solve a pre-image problem by getting back from the feature space to the input space. We also give new insights into the convexity of the pre-image problem. The relevance of the proposed method is evaluated on several time series.

INTRODUCTION

Many applications involve the analysis of time series data. One of the simplest, yet efficient, way to model time series is the autoregressive (AR) model. It states that each sample is given as a linear combination of a small number of previous samples. A prediction scheme is therefore inherent to this model, once the model parameters determined, i.e., the coefficients in the linear combination. These parameters are tightly linked to the covariance function of the process. The Yule-Walker equations explore this direct correspondence in order to estimate the parameters from the covariances of the time series. As a linear prediction model, the AR model is not adapted to treat nonlinear systems.

An elegant way to extend linear models into nonlinear ones is given by the concept of kernel methods in machine learning. The main idea is to map the data, using a nonlinear function, from the input space into a feature space, usually of a higher dimension. By using the kernel trick [START_REF] Aizerman | Theoretical foundations of the potential function method in pattern recognition learning[END_REF], it turns out that one can transform linear techniques into nonlinear ones, This work is partly supported by the Lebanese University and the French-Lebanese research program CEDRE No. 10 SCI F15/L5. without the need to explicitly exhibit the mapped space. This principle has shown its capacity in many applications, initially with Vapnik's Support Vector Machines (SVM) [START_REF] Vapnik | Statistical Learning Theory[END_REF], and now includes kernel principal component analysis and SVM novelty detection, only to name a few. In the same spirit, some kernel-based methods were considered for the analysis and prediction of time series data [START_REF] Richard | Online prediction of time series data with kernels[END_REF], including the SVM for regression and kernel Kalman filter [START_REF] Ralaivola | Time series filtering, smoothing and learning using the kernel kalman filter[END_REF].

In this paper, we explore the concept of kernel methods in order to provide a nonlinear extension of the AR model. To this end, we propose to take full advantage of the Yule-Walker equations in the feature space. We show that the parameters are estimated using the (lagged) expected kernels. The concept of expected kernels has shown its capacity in recent research [START_REF] Anderson | Channel-robust classifiers[END_REF][START_REF] Anderson | Expected kernel for missing features in support vector machines[END_REF]. Finally, to provide a prediction, one needs to map the result back into the input space. This is the pre-image problem. We give in this paper new insights into the convexity of this problem, and propose a technique to solve the problem. See [START_REF] Honeine | The pre-image problem in kernelbased machine learning[END_REF] for a recent review, with several applications in signal processing.

The rest of the paper is organized as follows: Next, we present the classical AR model with Yule-Walker equations. In Section 2, we present the kernel-based AR model and detail the solution of Yule-Walker equations in feature space. In Section 3, we provide a prediction scheme by solving the preimage problem. Section 4 covers the experiments done to evaluate the relevance of the proposed technique on "MG 30 " and "Lorenz attractor" time series.

AR model: the Yule-Walker equations

The AR model is widely used to analyze stationary and nonstationary time series [START_REF] Chevalier | Autoregressive models with time-dependent log area ratios[END_REF]. It gives each sample as a linear combination of previous samples. The AR model is defined by the fixed weights α i , for i = 1, 2, . . . , p, where p defines the order of the model. Let x 1 , x 2 , . . . , x n be a time series, a p-order AR model is described by

x i = p j=1 α j x i-j + ǫ i ,
for i = p + 1, . . . , n, where ǫ i is the unfitness error, assumed white Gaussian with zero mean. The parameters α 1 , α 2 , . . . , α p are directly connected with the covariance function of the process. One can therefore determine these parameters from the autocorrelation function. This is the essence of the Yule-Walker equations: Let r be the autocorrelation function of the time series, then r(τ ) = p j=1 α j r(τj), for a lag τ ≥ 1. Since r(-τ ) = r(τ ), we obtain the matrix form of the Yule Walker equations

r = R α, where r = [r(1) • • • r(p)] ⊤ , α = [α 1 • • • α p ] ⊤ , and R =      1 r(1) . . . r(p -1) r(1) 1 . . . r(p -2) . . . . . . . . . r(p -1) r(p -2) . . . 1     
, where r(0) = 1 without loss of generality. Assuming that the p × p matrix R is invertible, the coefficients α are estimated by α = R -1 r. Once the coefficients are estimated, the AR model can be applied to predict future samples.

KERNEL AUTOREGRESSIVE MODEL USING YULE-WALKER EQUATIONS

In order to derive a nonlinear extension of the Yule-Walker equations for autoregressive models, we use the principle of kernel machines. Let X be the input space, and let the kernel κ : X × X → IR be positive semi-definite, namely i,j α i α j κ(x i , x j ) ≥ 0 for all x i , x j ∈ X and any α i , α j ∈ IR. The Moore-Aronszajn theorem [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] states that a positive semi-definite kernel corresponds to an inner product in some arbitrary feature space. Let Φ(•) denotes the mapping function from the input space X into the feature space H, then κ(x i , x j ) = Φ(x i ), Φ(x j ) H , for any x i , x j ∈ X , where • , • H denotes the corresponding inner product in H, and • H its norm.

Using this concept, each sample of the time series x 1 , x 2 , . . . , x n is mapped from the input space into some feature space, yielding Φ(x 1 ), Φ(x 2 ), . . . , Φ(x n ). Thus, an AR model in the feature space is defined by

Φ(x i ) = p j=1 α j Φ(x i-j ) + ε Φ i ,
which belongs to the feature space span by the images of the samples with the map Φ(•). While the samples x i are assumed zero-mean, this is not often the case for the mapped data Φ(x i ). Let µ be the mean of the mapped time series, namely

µ = IE[Φ(x i )],
where IE[•] is the expectation, with i Φ(x i ). Then,

Φ(x i ) -µ = p j=1 α j Φ(x i-j ) + ε Φ i -µ = p j=1 α j Φ(x i-j ) -µ + ε Φ i -1 - p j=1 α j µ.
With the expectation of both sides, the times series being assumed stationary, we get (1

-p j=1 α j )µ = IE[ε Φ i ]
. By considering the inner product (in the feature space) of both sides of the above equation with (Φ(x i-τ )µ), for some positive lag τ , we get

Φ(x i ) -µ, Φ(x i-τ ) -µ H = ε Φ i -IE[ε Φ i ], Φ(x i-τ ) -µ H + p j=1 α j Φ(x i-j ) -µ, Φ(x i-τ ) -µ H . (1)
By substituting the inner product in the feature space with the corresponding kernel function, we define the centered version of a kernel κ(•, •) with

κ c (x i , x j ) = Φ(x i ) -µ, Φ(x j ) -µ H .
By analogy with the linear AR case, we assume that the noise ε Φ i and Φ(x i-τ ) are uncorrelated for every positive lag τ . Therefore, by taking the expectations of expression (1) and assuming the stationarity, we get for any τ ≥ 1:

IE[κ c (x i , x i-τ )] = p j=1 α j IE[κ c (x i-j , x i-τ )], (2)
where the notion of expected kernels is equivalent to the one recently studied in [START_REF] Anderson | Expected kernel for missing features in support vector machines[END_REF]. By considering all the lag values, expression (2) is written in matrix form

r κ = R κ α,
where

r κ = IE[κ c (x i , x i-1 )] IE[κ c (x i , x i-2 )] • • • IE[κ c (x i , x i-p )] ⊤ ,
and R κ is the matrix described by the expected kernels with

     IE[κ c (x i , x i )] IE[κ c (x i-2 , x i-1 )] • • • IE[κ c (x i-p , x i-1 )] IE[κ c (x i-1 , x i-2 )] IE[κ c (x i , x i )] • • • IE[κ c (x i-p , x i-2 )] . . . . . . . . . . . . IE[κ c (x i-1 , x i-p )] IE[κ c (x i-2 , x i-p )] • • • IE[κ c (x i , x i )]     
The vector of coefficients α is obtained by inverting the matrix R κ , with α = R -1 κ r κ .

In practice, the expectations are estimated over a set of n available samples. The centered version of the kernel is evaluated using

κ c (x i , x j ) = κ(x i , x j ) - 1 n n k=1 κ(x i , x k ) - 1 n n k=1 κ(x j , x k ) + 1 n 2 n k,k ′ =1 κ(x k ′ , x k ).

PREDICTION BY SOLVING THE PRE-IMAGE PROBLEM

Once the model parameters are determined on a set of n available samples, one may predict some future sample by

ψ i = p j=1 α j Φ(x i-j ),
starting with i = n + 1. Being a linear combination of p images, ψ i belongs to the feature space. In order to predict a sample, one needs to map back ψ i from the feature space into the input space. This is the pre-image problem. Several techniques have been proposed to solve this ill-posed problem. This includes a multidimensional scaling technique, a conformal map, and a learning scheme. See [START_REF] Honeine | The pre-image problem in kernelbased machine learning[END_REF] for a formal definition of the pre-image problem, and a recent survey of the literature.

We solve the pre-image problem by seeking an approximate solution x * i whose counterpart Φ(x * ) is as close as possible to ψ i , the latter being defined as above. The resulting optimization problem is:

x * i = arg min x 1 2 Φ(x) - p j=1 α j Φ(x i-j ) 2 H .
This is equivalent to the optimization problem

x * i = arg min x J i (x), with J i (x) = - p j=1 α j κ(x i-j , x) + 1 2 κ(x, x), (3)
where the term independent of x has been removed. Consider the case of the Radial Basis Functions, with kernels of the form

κ(x j , x j ′ ) = f ( x j -x j ′ 2 ). (4)
A sufficient condition for a function f ∈ C ∞ to be a valid positive-definite kernel is its complete monotonicity, i.e., its k-th derivative satisfies

(-1) k f (k) (ζ) ≥ 0 (5)
for any non-negative ζ [START_REF] Cucker | On the mathematical foundations of learning[END_REF]. The following Proposition gives insights on the convexity of the pre-image problem, as defined by (3):

Proposition 1. A sufficient condition for the convexity of the cost function is given by the non-negativity of the coefficients α 1 , . . . , α p .

Proof. Taking the second derivative of the cost function (3) with respect to x, we get

∇ 2 x J i (x) = ∇ x 2 p j=1 α j (x i-j -x) f (1) ( x i-j -x 2 ) = 2 p j=1 α j -f (1) ( x i-j -x 2 ) +2(x i-j -x) 2 f (2) ( x i-j -x 2 )
The term between parentheses is positive, due to condition [START_REF] Anderson | Channel-robust classifiers[END_REF]. Therefore, a sufficient condition for the second derivative to be positive, and thus for the convexity of (3), is that all the coefficients α j are positive.

Unfortunately, the cost function is not convex in the general case. However, it is reasonable to consider a local model, since the pre-image is estimated from its p previous samples. Therefore, we consider a local gradient descent approach with

x * i,t+1 = x * i,t -η t ∇ x J i (x * i,t ),
where the index t denotes the iterative technique. The convergence is controlled by the step size η t , in the opposite direction of the gradient of J i (x) with respect to x. Using the form (4), the gradient of the kernel with respect to x is given by

∇ x κ(x i-j , x) = 2(x i-j -x) f (1) ( x i-j -x 2 ).
By combining this expression with the gradient of the cost function J i (x), we get

∇ x J i (x) = 2 p j=1 α j (x i-j -x) f (1) ( x i-j -x 2 ).
Such expression simplifies further for several kernel functions, such as the Gaussian kernel with

f (ζ) = exp -1 2σ 2 ζ , thus f (1) (ζ) = -1 2σ 2 f (ζ).

EXPERIMENTS

We illustrate the efficiency of the proposed method with two well-known time series data: "Mackey-Glass" (MG 30 ) time series provides a model of the blood cells production evolution, and "Lorenz attractor" is the solution to a system of Fig. 1. Initial "Lorenz attractor" time series, with its prediction using the proposed approach. differential equations. We considered the Gaussian kernel. A set of n = 300 samples was used to determine the order p, the parameters α 1 , α 2 , . . . , α p of the AR model, as well as the value of the bandwidth σ of the Gaussian kernel. The next 300 samples, for i = n + 1, . . . , 2n, were used to evaluate the performance of the proposed kernel-based AR model. This configuration is identical to the one given in [START_REF] Ralaivola | Time series filtering, smoothing and learning using the kernel kalman filter[END_REF], from which we borrowed the comparative results given in Table 1. Depending on the bandwidth chosen for the Gaussian kernel, the error can be reduced to be approximately inconsiderable. As far as the bandwidth is getting smaller, the error will be decreased. The mean square error was estimated with

ǫ = 1 n 2n i=n+1 x * i -x i 2 ,
where x * i is the predicted value at instant i, and x i is the true value of the time series at the same time. An illustration of the prediction performance is given in Figure 1 for the "Lorenz attractor". It is worth noting that the proposed approach is simpler to implement, with less parameters, and significantly lower computational complexity than all the other methods given in Table 1.

CONCLUSION

We presented a kernel-based AR model for the prediction of time series data. We showed that one can take advantage of the Yule-Walker equations in the feature space, by using the concept of expected kernels. The prediction was assured by solving a pre-image problem. The relevance of this proposed method was revealed by a comparison to well-known nonlinear prediction techniques.

There are many possibilities for future work: we are currently working on methods to estimate the optimal order of the AR model in the feature space, in the same spirit of the Akaike Information Criterion. We are also studying other parameter estimation techniques, such as the Levinson-Durbin method. Also, we are considering the use of Gaussian process methods where the covariance function is the kernel function. Moreover, it would be desirable to extend the proposed approach to the autoregressive moving average model.

Table 1 .

 1 Mean square error for different nonlinear prediction approaches.

		"MG 30 "	"Lorenz"
	Multilayer perceptron	0.0461	0.2837
	Support vector regression	0.0313	0.1811
	Nonlinear Kalman filter	0.0307	0.03183
	Kernel AR model	0.00006	0.1793