Jie Chen
email: jie.chen@unice.fr

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

Paul Honeine
email: paul.honeine@utt.fr

Identification en ligne avec régularisation 1 Algorithme et analyse de convergence en environnement non-stationnaire

Cet article présente une méthode d'identification en ligne de système linéaire avec régularisation 1. L'analyse de convergence est effectuée pour des environnements non-stationnaires. Ce travail est une extension significative des algorithmes de la famille non-negative LMS, à la fois dans la forme de l'algorithme et de l'analyse de convergence. Par ses performances et son coût calculatoire réduit, l'algorithme présente des caractéristiques comparables à l'algorithme LMS tout en prenant en compte la régularisation 1. Des simulations valident la méthode proposée et le modèle de convergence.

Introduction

Several applications have recently shown the need for online sparse identification techniques. For instance, a particular driving force behind the development of such algorithms is the channel estimation problem, due to the fact that although the number of coefficients of the impulse is large, only a small portion has significant values. Compressive sensing theory provides a robust framework to estimate sparse signals. Instead of the accurate count of non-zero elements by 0 -norm, which leads to NP hard optimization problems, other sparse-induced norms can be used to overcome the difficulty caused by the 0 -norm. The use of the 1 -norm constraint or regularization is a popular choice.

Many approaches to solve the 1 -norm related problem have been described in the literature. Interior-point methods transfer these problems to a convex quadratic problem [START_REF] Nesterov | Interior-point polynomial methods in convex programming[END_REF]. Other recent methods include coordinate-wise descent methods [START_REF] Friedman | Pathwise coordinate optimization[END_REF], iterated shrinkage methods [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF], gradient methods [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF], gradient projection for sparse reconstruction algorithm [START_REF] Figueiredo | Gradient projection for sparse reconstruction : Application to compressed sensing and other inverse problems[END_REF], and Bregman iterative method [START_REF] Cai | Linearized bregman iterations for compressed sensing[END_REF].

However, the above methods all operate in batch mode. To identify sparse systems in an online way, several adaptive algorithms have also been proposed, including proportionate adaptive filters which incorporate the importance of the individual components by weights, 0 -norm constraint LMS algorithm which approximates 0 -norm by a differentiable function [START_REF] Gu | 0 -norm constraint lms algorithm for sparse system identification[END_REF], sparse LMS which uses the sign function as the subgradient of the 1 -norm [START_REF] Chen | Sparse lms for system identification[END_REF].

The last one of these three algorithms directly deals with system identification problems subject to 1 -norm constraint. However, the convergence of subgradient method is not guaranteed, even in batch mode, unless the step size is carefully chosen.

Considering that the 1 term can be rewritten as the sum of two non-negative vectors, the 1 -regularized problem can be transformed into a minimization problem under non-negativity constraints. Benefiting from this fact, we can solve this problem in an online manner efficiently using our previously proposed constrained system identification method, the so-called non-negative least-meansquare algorithm [START_REF] Chen | Non-negative least-mean-square algorithm[END_REF]. In this paper, we derive this algorithm and we propose models that characterize its meanweight behavior and its second-order behavior.

Presentation of the method

Consider an unknown linear system, parameterized by the discrete response vector α * of length N , with input x(n) and desired reference y(n). We intend to determine the coefficients of the system by minimizing the meansquare error with the sparsity-induced 1 -norm

α o = arg min α 1 2 E{[α x(n) -y(n)] 2 } + λ α 1 (1)
where the parameter λ provides a tradeoff between data fidelity and solution sparsity. This 1 -regularized problem can be easily rewritten as a standard non-negative leastsquare problem, by introducing two N -length non-negative vectors α + and α -which satisfy the following relations

α = α + -α - with α + , α -0 (2)
where is the component-wised symbol greater than or equal. These relations are satisfied by α + i = {α i } + and α - i = {-α i } -for all i = 1, 2, . . . , N, where {•} + denotes the positive-part operator defined as {x} + = max{0, x}. For simplicity, let us define a new vector α of length 2N by associating α + and α -as follows α = [α + α -] .

(3) We also define the extended input vector

x(n) = [x (n) -x (n)] . (4)
The problem (1) can thus be reformulated with respect to vector α by

αo = arg min α 1 2 E [ α x(n) -y(n)] 2 + λ 1 2N α subject to α 0 (5)
with 1 2N an all-one vector with 2N elements. Although the decomposition (2) is not unique, one can observe that the 2 -term is unaffected if we set α + ← α + + s and α -← α -+ s, where s 0 is a shift vector. However such a shift increases the regularization term in ( 5) by λ 1 2N s. It follows that, at the optimum for the problem (5), either

α + i = 0 or α - i = 0, for i = 1, 2, . . . , N so that in fact α + i = {α i } + and α - i = {-α i } + .
Note that the problem (5) has been reformulated as a system identification problem under non-negativity constraint with respect to α. Considering the stochastic gradient approximation in which the correlation matrix R x and cross-correlation r xy are replaced by their instantaneous estimates x(n)x (n) and x(n) y(n), we can now solve this problem in an online manner based on our previously proposed non-negative LMS algorithm [START_REF] Chen | Non-negative least-mean-square algorithm[END_REF] α

(n + 1) = (1 -η λ) α(n) + η D α(n) e(n) x(n) (6)
where η is the step size, D α(n) is the diagonal matrix with the ith diagonal element αi (n), and e(n) is the estimation error such that e

(n) = y(n) -α (n) x(n). Using the rela- tion (2), at each time instant the system coefficients α(n) is obtained by α(n) = α + (n) -α -(n) (7)
This algorithm is based on the NNLMS algorithm. Thus, the variants of NNLMS can also be applied to improve its performance [START_REF] Chen | A modified non-negative lms algorithm and its stochastic behavior analysis[END_REF].

Algorithm behavior modeling in non-stationary environments

We now study the mean-weight behavior of the proposed adaptive algorithm [START_REF] Figueiredo | Gradient projection for sparse reconstruction : Application to compressed sensing and other inverse problems[END_REF] in an arbitrary but quite general time-variant environment. The input x(n) and the desired output y(n) signals are assumed stationary and zeromean. The signal z(n) = y(n)x α * accounts for measurement noise and modeling errors. It assumed that z(n) is stationary, zero-mean with the variance σ 2 z and statistically independent of any other signal. The time-variant environment is defined with respect to the system coefficients α * as α

* (n) = α * o (n) + ξ(n) (8)
where α * o (n) is a deterministic time-variant mean, and ξ(n) is a zero-mean random variable with covariance Ξ = σ 2 ξ I and independent of any other signal. This simple model provides some information on how the performance of the proposed algorithms is affected by a time-variant optimal solution which consists of a deterministic trajectory and a random perturbation. The model ( 8) leads to a tractable analysis and permits inferences about the behavior of the algorithms in time-variant environments by varying the mean value α * o (n) and the power σ 2 ξ of ξ(n).

Mean weight behavior analysis

We define the weight-error vector with respect to the mean coefficients α * o (n) by

ṽ(n) = α(n) -α * o (n) (9)
For the feasibility of the analysis, the following independence assumption is considered in the derivation.

Assumption 1 (Independence assumption) The input signal x(n) is independent of the weight error vector ṽ(m), for all time index m ≤ n.

Although not true in general, this assumption is commonly used for adaptive filter analysis, and the analytical result is usually not sensitive to this approximation. Now using the relation [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] in the update equation ( 6), we get an update equation for the weight-error vector

ṽ(n + 1) = (1 -η λ)ṽ(n) -∆(n) + η D x(n) α(n)e(n) (10)
where

∆(n) = α * o (n + 1) -(1 -η λ) α * o (n).
Note that the estimation error can be expressed by

e(n) = z(n) -ṽ (n)x(n) + ξ (n)x(n)
Now taking the expectation of [START_REF] Nesterov | Interior-point polynomial methods in convex programming[END_REF], neglecting the statistical dependence of x(n) and ṽ(n), observing that the vector ξ(n) is zero-mean and independent of the other signals, and using E {z(n)D x(n)} = 0, yields the meanweight behavior model

E {ṽ(n + 1)} = (1 -η λ)I -η D α * o (n) Rx E {ṽ(n)} -∆(n) -η diag{ Rx K(n)} (11)
where K(n) = E{ṽ(n)ṽ (n)} is the covariance matrix of ṽ(n), and Rx is the covariance matrix of x(n). Note that the recursion (11) needs the second-order statistics of ṽ(n). In order to simplify the model, we can make the following assumption K(n) ≈ E{ṽ(n)}E{ṽ (n)} in (11) to provide a simplified model. The first order analysis of adaptive filters is usually insensitive to this kind of approximation. As to be seen in experiments, the simplified model will provide accurate behavior description.

Excess mean-square error analysis

Neglecting the statistical dependence of x(n) and ṽ(n), and using the properties assumed for z(n) and ξ(n), yields an expression for the mean-square error model

E e 2 (n) = E (z(n) -ṽ (n) x(n) + ξ (n) x(n)) 2 = σ 2 z + trace{ Rx K(n)} + trace{R x Ξ}
The excess mean-square error (EMSE), which is usually more favorable in the performance analysis, is correspondingly provided by

J EMSE (n) = trace{ Rx K(n)} + trace{R x Ξ} (12)
The term trace{R x Ξ} is the contribution of the nonstationarity of the system to the excess error caused by random perturbation. In order to determine the excess mean square error due to trace{ Rx K(n)}, we determine a recursion for K(n), which is provided in Appendix due to its complexity. This result can be used to study the convergence behavior of E{e 2 (n)} and J EMSE (n).

Experiments

In order to validate the proposed algorithm and the models, we considered non-stationary environments defined by adding two types of time variant terms to stationary coefficients, denoted by "nst.1" and "nst.2", respectively. The system coefficients in these two environments are given by α * (nst.1) i

(n) = α * (st.) i

+ ξ i (n) and α * (nst.2) i (n) = α * (st.) i + |α * (st.) oi | 10 sin 2 π T n+2 π i -1 N +ξ i (n)
where α * (st.) in the above two expressions is the stationary component defined by

α * (stat.) i =      0.55 -0.1 i i= 1, . . . , 5 0 i = 6, . . . , 25 -0.1 (i -25) i = 26, . . . , 30
In this first case, there is only one random perturbation added to the stationary unconstrained solutions. Whereas a deterministic sinusoidal time-varying trajectory is also added in the second case. The period T of sinusoidal components was set to 2500. The input signal was an AR process given by x(n) = 0.5 x(n -1) + w(n), with w(n) i.i.d. zero-mean Gaussian with variance σ 2 w , adjusted to obtain the desired input power σ 2

x = 1. The step size was set to η = 0.005 and regularization parameter was set to λ = 0.06. The variance of modeling noise error z(n) remained σ 2 z = 0.01. Experiment results are shown in Fig. 1. The blue and red curves show the simulation results, and the theoretical predictions by (11) for the mean weight behavior, and by (12) for the EMSE . Simulation results were obtained by averaging 100 rMonte-Carlo runs. The mean weight behavior curves are illustrated in Figs. 1(a) and 1(b).

For the second-order EMSE curves, in addition to these curves, the variance σ 2 ξ was varied using the values in {0, 0.001, 0.005}. Corresponding EMSE curves obtained from the theoretical model are illustrated in Figs. 1(c) and 1(d). The simulation results conform with these curves but are not shown for clarity. Effects of the deterministic time-varying trajectory and random perturbations can be clearly observed in these figures. Extra EMSE arises due to tracking of the optimal solution variations.

Conclusion

In this paper, we studied the online system identification problem regularized by 1 -norm. This was performed by extending the NN-LMS algorithm. Analytical behavior model of the proposed algorithm was studied in nonstationary environments. Future work may include exploring variants of the algorithm to improve its performance.

Appendix : Recursion of K(n)

Post-multiplying [START_REF] Nesterov | Interior-point polynomial methods in convex programming[END_REF] by its transpose, taking the expected value, and using the statistical properties of z(n) and ξ(n), yields

K(n + 1) = (1 -ηλ) 2 K(n) -(P 1 (n) + P 1 (n)) -η (P 2 (n) + P 2 (n)) -η (P 3 (n) + P 3 (n)) + P 4 (n) + (P 5 (n) + P 5 (n)) + (P 6 (n) + P 6 (n)) + η 2 P 7 (n) + η (P 8 (n) + P 8 (n)) + η 2 P 9 (n) +
where matrices P 1 (n) to P 12 (n) are calculated by

P 1 (n) = (1 -ηλ)E {v(n)} ∆ (n) P 2 (n) = (1 -ηλ) K(n) Rx D ṽ (n) P 3 (n) = (1 -ηλ) K(n) Rx D α * o (n) P 4 (n) = E{∆(n) ∆ (n)} P 5 (n) = ∆(n) diag{ Rx K(n)} P 6 (n) = ∆ E{ṽ (n)} Rx D α * o (n) P 7 (n) = σ 2 z ( Rx • K(n)) P 8 (n) = σ 2 z E {D ṽ (n)} Rx D α * o (n) P 9 (n) = σ 2 z D α * o (n) Rx D α * o (n) P 10 (n) = Q(n) • K(n) P 11 (n) = E {D ṽ (n)} Q(n) D α * o (n) P 12 (n) = D α * o (n) Q(n) D α * o (n) where in the above expressions Q(n) = trace{ Rx K(n)} Rx + 2 Rx K(n) Rx .
These terms do not involve the random perturbation ξ(n) on the weights. We now derive expressions for P 13 (n) through P 15 (n). These terms convey the effect of the environment non-stationarity due to the random variations of system weights. These terms are expressed by

P 13 (n) = σ 2 ξ K(n) • ( Rx trace{R x } + 2 R x R x ) P 14 (n) = σ 2 ξ α * o (n) α * o (n) •( Rx trace{R x } + 2 R x R x ) P 15 (n) = σ 2 ξ E {ṽ(n)}α * o (n) •( Rx trace{R x }+2 R x R x ) with R x = [R x , -R x ] .
Using the expected values P 1 (n) to P 15 (n) in (13), we finally obtain a recursive analytical model for the behavior of K(n). More detailed derivations are similar to those of NN-LMS algorithm and its variants, and can be referred to [START_REF] Chen | A modified non-negative lms algorithm and its stochastic behavior analysis[END_REF][START_REF] Chen | Non-negative least-mean-square algorithm[END_REF].

  2nd-order behavior with α * (nst.1) 2nd-order behavior with α * (nst.2)

Figure 1 -

 1 Figure 1 -Algorithm behavior in non-stationary environments. First row : first order results with σ 2 ξ = 0.005. Second row : Second order results with various levels of σ 2 ξ .


N Ghadban

P Honeine
email: paul.honeine@utt.fr

C Francis
email: cfrancis@ul.edu.lb

F Mourad-Chehade

J Farah

M Kallas
email: maya.kallas@utt.fr

Estimation locale d'un champ de diffusion par modèles à noyaux

Introduction

Les Réseaux de Capteurs sans Fil (RCF) sont introduits comme une solution économique et facilement déployable pour la surveillance à distance de phénomènes physiques [START_REF] Akyildiz | A survey on sensor networks[END_REF]. Initialement développés pour des applications militaires, leur domaine d'application ne cesse de s'élargir pour inclure la surveillance dans les domaines civil, industriel et environnemental [START_REF] Steere | Research challenges in environmental observation and forecasting systems[END_REF]. Les RCF sont constitués de noeuds déployés en grand nombre en vue de collecter et de transmettre les données environnementales, d'une manière autonome. Cette autonomie impose des contraintes d'énergie et de débit de communication, ce qui exige une optimisation des algorithmes correspondants : les capteurs communiquent de proche en proche, entre voisins. Dans ce papier, nous traitons le problème de la modélisation de la diffusion d'un champ de gaz, ainsi que le suivi de son évolution au cours du temps. Notre but revient ainsi à estimer la distribution de gaz aux endroits dépourvus de capteurs. Vu l'application considérée, nous devons recourir à des méthodes de modélisation adéquates. Pour cela, nous considérons un cadre similaire à celui donné dans [START_REF] Honeine | A decentralized approach for non-linear prediction of time series data in sensor networks[END_REF], tout en proposant une approche locale pour l'estimation.

Description de la méthode

Supposons que le réseau est composé d'un grand nombre de capteurs déployés dans une région χ, où χ ⊂ R 2 pour un espace de dimension 2, ou bien χ ⊂ R 3 comme dans le cas des drones navigant dans un espace de dimension 3. Les capteurs communiquent entre eux de proche en proche. Considérons un capteur i de position x i ∈ χ, et soit ϑ i l'ensemble des indices des k plus proches capteurs voisins au capteur i. Dans notre approche, nous considérons que x i est voisin à lui-même, c'est à dire i ∈ ϑ i . A chaque instant, chaque capteur i reçoit de ses voisins j ∈ ϑ i , j = i, des informations comportant leurs positions x j , et les valeurs de la quantité physique mesurée y j . L'objectif est alors de trouver, pour chaque capteur, un modèle local ψ i (•) pour représenter le champ mesuré de façon à avoir ψ i (x j ) = y j pour tout j ∈ ϑ i . Nous proposons d'utiliser le formalisme des méthodes à noyaux [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] pour trouver ces modèles.

Le problème de minimisation de l'erreur quadratique entre les sorties du modèle ψ i (x j ) et les réponses mesurées y j nous amène au problème d'optimisation suivant :

ψ * i = arg min ψ∈H j∈ϑi (y j -ψ i (x j )) 2 + η ψ i 2 H ,
où H est un espace de Hilbert à noyau reproduisant κ(•, •) constitué de fonctions de χ dans R, η est un paramètre de régularisation qui contrôle le compromis entre le raccord aux données 1 disponibles et la douceur de la solution. D'après le théorème de représentation [START_REF] Schölkopf | A Generalized Representer Theorem[END_REF], la solution de ce problème est donnée sous la forme d'une combinaison linéaire de noyaux κ(x j , •), selon

ψ * i (•) = j∈ϑi α i,j κ(x j , •),
où les coefficients α i,j , pour j ∈ ϑ i , sont à déterminer. Soient α i le vecteur contenant ces coefficients et y i le vecteur des mesures effectuées y j , pour j ∈ ϑ i . Nous démontrons que :

α i = (K i + ηI) -1 y i ,
où K i est la matrice d'éléments κ(x m , x n ), pour m, n ∈ ϑ i , et I la matrice identité de même taille.

Simulations

Afin d'illustrer les résultats de ce travail, nous considérons la diffusion d'un gaz dans un espace carré χ = [-0.8; 0.8] × [-0.8; 0.8]. Une source de gaz placée au point S(0.47; 0.42) est ainsi activée de t = 1 jusqu'à t = 100. Les capteurs déployés aléatoirement fournissent des mesures relevées en tout instant t. Sur 100 capteurs, 66 sont utilisés pour l'apprentissage et 34 pour la validation du modèle calculé. Soit Ω val l'ensemble de validation. La Figure 1 représente la répartition des capteurs dans χ. Les cercles bleus désignent les capteurs d'apprentissage, et les carrés rouges, numérotés de 1 à 34, désignent les capteurs de validation. Nous utilisons dans ce papier le noyau Gaussien donné par κ(x m , x n ) = exp(-1 2σ 2 x mx n 2 ), σ étant la largeur de bande du noyau. La Figure 2 illustre les résultats obtenus. On remarque un bon accord entre les mesures du champ de gaz réel et celles estimées à partir du modèle.

FIG. 1 -

 1 FIG. 1 -Répartition des capteurs (apprentissage • et validation ).

FIG. 2 -

 2 FIG.2 -a-Estimation au capteur 3 ; b-Erreur quadratique mesurée sur l'ensemble Ω val ; c-Erreur absolue au capteur 1.


Jie Chen

André Ferrari
email: andre.ferrari@unice.fr

Paul Honeine
email: paul.honeine@utt.fr

NONLINEAR UNMIXING OF HYPERSPECTRAL DATA WITH PARTIALLY LINEAR LEAST-SQUARES SUPPORT VECTOR REGRESSION

Keywords: Nonlinear unmixing, hyperspectral image, support vector regression, multi-kernel learning, spatial regularization

In recent years, nonlinear unmixing of hyperspectral data has become an attractive topic in hyperspectral image analysis, because nonlinear models appear as more appropriate to represent photon interactions in real scenes. For this challenging problem, nonlinear methods operating in reproducing kernel Hilbert spaces have shown particular advantages. In this paper, we derive an efficient nonlinear unmixing algorithm based on a recently proposed linear mixture/nonlinear fluctuation model. A multi-kernel learning support vector regressor is established to determine material abundances and nonlinear fluctuations. Moreover, a low complexity locally-spatial regularizer is incorporated to enhance the unmixing performance. Experiments with synthetic and real data illustrate the effectiveness of the proposed method.

INTRODUCTION

Hyperspectral imagery provides data with a wide spectral range and high spectral resolution. These characteristics are desirable for detection and classification of surfaces and chemical elements in real images. Hyperspectral image processing has become a continuously growing area, and now receives considerable attention.

Due to multiple factors, including the possible low spatial resolution of some hyperspectral-imaging devices, the diversity of materials, the reflections of photons onto several objects, etc., problems with mixed pixels can occur and be critical for proper interpretation of images. Spectral unmixing is an important issue to analyze hyperspectral data. It consists of decomposing each mixed pixel into pure endmember spectra, and estimating the fractional abundance of each endmember [START_REF] Keshava | Spectral unmixing[END_REF].

Linear unmixing methods have been extensively considered within the area of remote sensing. In the recent years, researchers have progressively focused on nonlinear unmixing method, as there are many situations in which the classical linear mixture model may be insufficient [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: geometrical, statistical, and sparse regression-based approches[END_REF]. Recent nonlinear unmixing approaches dedicated to hyperspectral images can be divided into the following classes: 1) Direct methods: Some techniques mimic the nonlinear interactions between endmember spectra by incorporating them into the regression model in an explicit way. For instance, the authors 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] extended the collection of endmembers by adding artificial crossed spectra to model photon interactions. 2) Manifold learning methods: These approaches exploit low dimensional manifolds embedded in high-dimensional sensed data [START_REF] Bachmann | Exploiting manifold geometry in hyperspectral imagery[END_REF]. Algorithms based on Isomap and LLE, which assume the local linearity of data, were recently used. 3) Kernel-based methods: These approaches can model complex interactions among material signatures by the virtue of kernel trick.

In [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinearfluctuation model[END_REF], we introduced a new framework for estimating abundances involved in a nonlinear mixture of hyperspectral data. It overcomes the drawbacks of former heuristic kernel unmixing methods [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF]. 4) Bayesian methods: In [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF][START_REF] Dobigeon | Supervised nonlinear spectral unmixing using a polynomial post nonlinear model for hyperspectral imagery[END_REF], Bayesian inference was used to estimate the fraction of abundances within the context of nonlinear mixing scenarios, using generated samples distributed according to the posterior distribution. 5) Neural networks: This family refers to techniques that intend to establish, based on a collection of training data, an input-output relation between some observed data and their abundances via black-box modeling [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].

Among the above mentioned methods, kernel-based algorithms have demonstrated great generality and excellent performance for nonlinear unmixing [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinearfluctuation model[END_REF][START_REF] Chen | A novel kernel-based nonlinear unmixing scheme of hyperspectral images[END_REF][START_REF] Chen | Nonlinear unmixing of hyperspectral images with multi-kernel learning[END_REF]. In this paper, the inherent spectral mixing mechanism is modeled by a partial linear model, consisting of linear mixture components and a nonlinear fluctuation term defined in a reproducing kernel Hilbert space (RKHS). A multi-kernel support vector regression algorithm is derived to estimate the material abundances, and to automatically adjust the proportion between the linear and the nonlinear components. Furthermore, spatial correlation is also considered by a local regularization term which does not introduce extra computational complexity. This work is a followon work of the state-of-the-art kernel unmixing work [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinearfluctuation model[END_REF][START_REF] Chen | Nonlinear unmixing of hyperspectral images with multi-kernel learning[END_REF]. Benefiting from the specific structure of the unmixing problem, and from the spatial information, the proposed approach provides a reliable estimation of fractional abundances. Its advantages are illustrated via comparative experiments.

NOTATIONS AND NONLINEAR MIXTURE MODEL

Consider that an observed scene consists of R significant endmembers with spectral signature mi ∈ R L , where L denotes the number of spectral bands. Note that, usually, we have R ≪ L. Let r ∈ R L be an hyperspectral pixel, and α ∈ R R the unknown abundance vector associated to the latter. Let M = [m1, . . . , mR] ∈ R L×R be the matrix of the endmember spectra. For the sake of convenience, the ℓ-th row of M is denoted by m ⊤ λ ℓ ∈ R L , that is, m λ ℓ is the vector of the endmember signatures at the ℓ-th wavelength band.

Nonlinear mixture usually results from complex interactions of light sources scattered by multiple materials in the scene. Linear mixing model is a reasonable approximation in cases where these interactions are not significant. However, when such effects cannot be neglected, it is more reasonable to assume that mixing consists of a linear trend and a nonlinear fluctuation term, rather than leaving a nonlinear model totally uncontrolled. Based on this point of view, we model the mixing mechanism with respect to each spectral band

ψ(m λ ℓ ) = θ α ⊤ m λ ℓ + ψnlin(m λ ℓ ) + e ℓ ( 1
)
where θ is an additional factor that affects the linear mixture, which is parameterized by the abundance vector α. Function ψnlin is supposed to describe the nonlinear interactions between spectra, and e ℓ is the modeling error. In this paper, the function ψnlin is restricted to be an element of a RKHS, say Hnlin, such that

ψnlin(m λ ℓ ) = ψ, κnlin(•, m λ ℓ ) H nlin , ∀ψnlin ∈ Hnlin (2)
The corresponding kernel satisfies to the reproducing property

κnlin(m λ ℓ , m λ k ) = κnlin(•, m λ ℓ ), κnlin(•, m λ k ) (3)
Several nonlinear models, for instance the bilinear model and postnonlinear models, can be expressed or approximated by this model. Moreover, the non-negativity and the sum-to-one constraints are usually considered for abundance α estimation, namely, α 0 and

1 ⊤ α = 1 (4)
It should be noticed that the rationality of imposing the sum-to-one constraint is under discussion among researchers. In this paper, we preserve it because discarding this constraint will only simplify the mathematical problem and the proposed algorithm. They can be easily modified to be adapted to this case.

NONLINEAR UNMIXING WITH LEAST-SQUARES SUPPORT VECTOR REGRESSION

We assume that the endmember matrix M has been determined by an endmember extraction algorithm. This approach is usually referred to as supervised unmixing. Pixel unmixing using the model ( 1) is achieved by estimating θ, α, and ψnlin. In order to further control the balance between the linear trend and the nonlinear fluctuation term, we propose to solve the following multi-kernel least-squares support vector problem, in which these proportions are controlled via functional norms.

ψ * ,θ * ,u * = arg min ψ,θ,u 1 2 ψlin 2 H ′ lin u + ψnlin 2 H ′ nlin 1 -u + 1 2µ L ℓ=1 e 2 ℓ
subject to e ℓ = r ℓψ(m λ ℓ )

and ψ = ψlin + ψnlin with ψlin(m λ ℓ ) = θ α ⊤ m λ ℓ α 0 1 ⊤ α = 1 θ ∈ IR * + 0 ≤ u ≤ 1
where Hlin and Hnlin are RKHS. This problem is not convex, and is difficult to solve as formulated. Fortunately, it can be transformed into an equivalent convex problem in the sense of [14, p. 130]. Consider the variable change h = θ α. The cost function can be directly reformulated as a function of h. The linear function is expressed by

ψlin(m λ ℓ ) = h ⊤ m λ ℓ , (5)
and the two constraints over α become

(h 0) and (1 ⊤ h = θ with θ ∈ IR * + )
The latter constraint can be eliminated because it is trivial due to the former one. Because h = θ α and 1 ⊤ α = 1, we have α * = h * /θ * with θ * = 1 ⊤ h * . This leads us to the following problem

min u J(u) subject to 0 ≤ u ≤ 1 (6)
where J(u) is defined by

J(u) =              min ψ F (u, ψ) = 1 2 ψ lin 2 H lin u + ψ nlin 2 H nlin 1-u + 1 2µ L ℓ=1 e 2 ℓ subject to e ℓ = r ℓ -ψ(m λ ℓ ) ψ = ψlin + ψnlin ψlin(m λ ℓ ) = h ⊤ m λ ℓ h 0 (7)
Appropriate definitions for Hlin and Hnlin, and continuity considerations, are used for defining J(u) in u = 0 and 1. See [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinearfluctuation model[END_REF] for details. It can be shown that the problem ( 6)-( 7) is a convex optimization problem by virtue of the convexity of the so-called perspective function defined by f (u, ψ) = ψ 2 H lin/nlin /u over IR+ × Hlin/nlin [START_REF] Boyd | Convex Optimization[END_REF]. This allows to formulate a two-stage optimization procedure, with respect to ψ and u successively.

Solving with respect to ψ

By the strong duality property, we shall now derive a dual problem that has the same solution J(u) = F (u, ψ * ) as the primal problem [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF]. Let us introduce the Lagrange multipliers β ℓ and γr. The Lagrange function associated with the problem [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF] can be written as

G = 1 2 1 u h 2 + 1 1 -u ψnlin 2 H nlin + 1 2µ L ℓ=1 e 2 ℓ - L ℓ=1 β ℓ (e ℓ -r ℓ + ψ(m λ ℓ )) - R r=1 γr hr (8)
with γr ≥ 0, where we have used that ψlin 2 H lin = 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)
Substituting ( 9) into (8) leads us to the dual problem [START_REF] Guilfoyle | A quantitative and comparative analysis of linear and nonlinear spectral mixture models using radial basis function neural networks[END_REF], see next page, where Knlin is the Gram matrix defined by

[Knlin] ℓk = κnlin(m λ ℓ , m λ k ) (11)
Problem ( 10) is a standard quadratic programming problem which can be solved by a generic solver. Once β * and γ * are determined, pixel reconstruction can be performed using [START_REF] Dobigeon | Supervised nonlinear spectral unmixing using a polynomial post nonlinear model for hyperspectral imagery[END_REF]. Finally, using the relation between h and α, the estimated abundance vector is given by

r * = [ψ * (m λ 1 ), . . . , ψ * (m λ L )] ⊤ with ψ * (m λ ℓ ) = m ⊤ λ ℓ h * + ψ * nlin (m λ ℓ ) defined in equation
α * = M ⊤ β * + γ * 1 ⊤ (M ⊤ β * + γ * ) (12)

Solving with respect to u

In the general framework for multi-kernel learning described in [START_REF] Rakotomamonjy | Sim-pleMKL[END_REF], the kernel weights were updated by the reduced gradient method.

In [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinearfluctuation model[END_REF][START_REF] Chen | Nonlinear unmixing of hyperspectral images with multi-kernel learning[END_REF], a gradient descent was also used to update u. However, gradient based updates usually have low convergence rates. An optimal step size search strategy is often considered to accelerate the convergence rate, with a an additional computational cost. We shall now show that there is a closed-form expression for the optimum u * . Notice that the function fp,q(u) = p u + q 1u with p, q ≥ 0 (13)

J(u) =          max β,γ G ′ (u, β, γ) = - 1 2 β γ ⊤ Ku + µI uM uM ⊤ uI β γ + r 0 ⊤ β γ subject to γ 0 with Ku = uM M ⊤ + (1 -u)Knlin (10)
Jsp(u) =              max β,γ G ′ sp (u, β, γ) = - 1 2 β γ ⊤ Ku + µI ξM ξM ⊤ ξI β γ + r -ξζ M i∈Nn ωih * i -ξζ i∈Nn ωih * i ⊤ β γ subject to γ 0 with Ku = ξM M ⊤ + (1 -u)Knlin and ξ = u 1 + uζ (19)
is convex over the interval (0, 1). The optimum value is achieved at

u * = (1 + q/p) -1 (14)
Now considering the problem (6), and using the stationarity conditions ( 9), we conclude that the optimum value u * at each iteration is given by

u * = 1 + (1 -u * -1 ) (β * ⊤ Knlinβ * )/ h * 2 -1 (15)
where u * -1 denotes the optimal u * obtained at the previous iteration.

LOCALLY-SPATIAL REGULARIZATION

One of the distinguishing properties of remotely sensed data is that they convey multivariate information into a 2D pictorial representation [START_REF] Plaza | Foreword to the special issue on spectral unmixing of remotely sensed data[END_REF]. Incorporating spatial information may have a positive impact on hyperspectral unmixing. Exploiting spatial information can be conducted by designing appropriate criteria, usually, by adding a penalizing term to the modeling error that promotes similarity of fractional abundances between one pixel and some others to be defined. Spatial regularization techniques incorporating all the pixels, such as TV norm, lead to a global optimization problem that notably increase the computational burden [START_REF] Iordache | Total variation spatial regularization for sparse hyperspectral unmixing[END_REF][START_REF] Chen | Nonlinear estimation of material abundances of hyperspectral images with l1-norm spatial regularization[END_REF]. To overcome this drawback, we adopt hereafter a local regularization strategy, and we associate it to the nonlinear unmixing algorithm presented before. Let us identify the current pixel under analysis by the index n. In addition to estimating hn via the optimization problem ( 6)-( 7), we shall also promote the similarity between hn and its neighbors by considering the quadratic penalty term

Jreg(hn) = ζ 2 i∈Nn ωi hn -h * i 2 . ( 16
)
where the dissimilarities with the neighboring pixels i in Nn are weighted by ωi. The abundance vectors h * i are considered to have already been estimated. If the image is treated one pixel after another, sequentially along rows, an example of such a neighborhood is defined by {n -1, nw, nw -1}, where w is the width of the image. Although ℓ1-type regularization has favorable properties in image processing, it adds significant complexity to the nonlinear algorithm. For this reason, we shall employ this ℓ2 regularization combined with the following strategy 1. Calculate ωi = rn -ri 2 / rn 2 , which represents the normalized spectral distance between pixel n and its neighbors;

2. If mini ωi > ν0, with ν0 a given threshold, the algorithm is conducted without spatial regularization using [START_REF] Guilfoyle | A quantitative and comparative analysis of linear and nonlinear spectral mixture models using radial basis function neural networks[END_REF].

3. Otherwise, the weights ωi are determined from the ωi's by

ωi = 1/ ωi i∈Nn 1/ ωi (17)
and the algorithm is conducted with the spatial regularization term (16).

Step 2. allows to disable spatial regularization in the case where the spectral distance between pixel n and its neighbors is significant, which may preserve transient phenomena. This step can obviously be eliminated if preferred, and the definition (17) changed. Finally, unmixing with spatial regularization is performed by adding ( 16) to the Lagrange function [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF]. The conditions for optimality yield

     h * n = u 1+uζ L ℓ=1 β * ℓ m λ ℓ + γ * + ζ i∈Nn ωih * i ψ * nlin = (1 -u) L ℓ=1 β * ℓ κnlin(•, m λ ℓ ) e * ℓ = µ β * ℓ (18)
Substituting [START_REF] Chen | Nonlinear estimation of material abundances of hyperspectral images with l1-norm spatial regularization[END_REF] into the Lagrange function leads us to the dual problem [START_REF] Baldridge | The ASTER spectral library version 2.0[END_REF]. It is interesting to note that this optimization problem has the same computational complexity as [START_REF] Guilfoyle | A quantitative and comparative analysis of linear and nonlinear spectral mixture models using radial basis function neural networks[END_REF]. Once the duality variables are determined, the vector of fractional abundances can be estimated as follows

α * n = M ⊤ β + γ * + ζ i∈Nn ωih * i 1 ⊤ M ⊤ β + γ * + ζ i∈Nn ωih * i . (20)
Finally, note that the update equation ( 15) for u remains unchanged.

EXPERIMENT RESULTS

Experiments with synthetic data

Let us first report some experimental results on synthetic images. The endmembers were randomly selected from the library ASTER where signatures have reflectance values measured over 224 spectral bands, uniformly distributed in the interval 3 -12 µm [START_REF] Baldridge | The ASTER spectral library version 2.0[END_REF].

Spatially-uncorrelated pixels

Random pixels were first generated to test the performance of the basic nonlinear unmixing algorithm. Pixels with three, five and eight endmembers were generated respectively with the linear model, the generalized bilinear mixture model [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF], and the post-nonlinear mixing model (PNMM) defined by r = (M α) 0.7 +n. In each scenario, 1000 pixels were generated for evaluating and comparing the performance. All these images were corrupted with an additive white Gaussian noise n with SNR = 30 dB. The Gaussian kernel defined as κnlin(m λp , m λ ℓ ) = exp(-m λpm λ ℓ 2 /2σ 2 0 ) was used. 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], and the extended endmember-matrix method (ExtM) [START_REF] Raksuntorn | Nonlinear spectral mixture analysis for hyperspectral imagery in an unknown environment[END_REF], and the K-Hype algorithm [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinearfluctuation model[END_REF]. The root mean square error (RMSE) of the abundance vectors was used to compare these algorithms. Preliminary runs were performed with independent data to set their parameters. The results are reported in Table 1. The proposed method had better performance than the other algorithms except for linearly-mixed data. It however outperformed K-Hype in this case because it can automatically balance ψlin and ψnlin. This is a key advantage because nonlinear mixing effects can be weak in some situations, in particular if there is a locally-dominant material. Figure 1 illustrates this situation, where material #9 is dominant in the center of the scene. See image DC2 in [START_REF] Iordache | Total variation spatial regularization for sparse hyperspectral unmixing[END_REF], mixed with a bilinear model.

Spatially-correlated pixels

We shall now check the advantage of incorporating spatial information, and validate the proposed locally-spatial regularization. A synthetic image of 75 × 75 pixels was generated with the bilinear mixture model using five signatures. The pure regions and mixed regions involved between 2 and 5 endmembers, spatially distributed in the form of square regions. See image DC1 in [START_REF] Iordache | Total variation spatial regularization for sparse hyperspectral unmixing[END_REF] for details. The image was corrupted by a white Gaussian noise with SNR = 25 dB. The regularization parameter ζ was set to 10. The threshold ν0 was set to 0.01. The true abundances and estimated abundances of the 4-th and 5-th materials are illustrated in Figure 2, with RMSE at the bottom of each image. The effect of spatial regularization can be clearly observed.

Experiments with real data

This experiment illustrates the performance of the proposed algorithms, when applied to real hyperspectral data. This scene is the well-known image captured on the Cuprite mining district (NV, USA) by AVIRIS. A sub-image of 250 × 191 pixels was chosen to evaluate the algorithms. This area of interest has L = 188 spectral bands. The number of endmembers was first estimated via the virtual dimensionality, and R was accordingly set to 12 [START_REF] Nascimento | Vertex Component Analysis: A fast algorithm to unmix hyperspectral data[END_REF]. VCA algorithm was then used to extract the endmembers. The proposed algorithms, with/without spatial regularization, were run with the Gaussian kernel. The kernel bandwidth was set to σ 2 0 = 2. The regularization parameter µ was fixed to 5 • 10 -3 . The spatial regularization parameter ζ was set to 10. The unmixing results for three materials are shown in Figure 3. It can be observed that the mineral distributions are consistent with the reference map [START_REF] Clark | Imaging spectroscopy: Earth and planetary remote sensing with the usgs tetracorder and expert systems[END_REF], and the US Highway 95 can be clearly identified whereas the linear unmixing failed to separate it [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinearfluctuation model[END_REF].

CONCLUSION AND PERSPECTIVES

We presented a least-squares support-vector regression approach for hyperspectral data unmixing. The model is partially linear, that is, it combines a linear trend and a nonlinear correction term. Multikernel learning was used to tune the balance between the linear and nonlinear components, resulting in a convex optimization problem.

A locally-spatial regularizer, which does not introduce extra computational expenses, was proposed to enhance the unmixing performance. Future work will focus on endmember extraction within this nonlinear unmixing framework.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Illustration of the situation where one material dominates the others. True Proposed Prop. + Reg.

Fig. 3 .

 3 Fig. 3. Unmixing results for the Cuprite district scene

Table 1 .

 1 Spatially-uncorrelated pixels: RMSE comparison (×10 -2 )

			R = 3			R = 5			R = 8
		linear	bilinear	PNMM	linear	bilinear	PNMM	linear	bilinear	PNMM
	FCLS	1.10 ± 0.02	33.2 ± 4.37	21.3 ± 1.31	1.65 ± 0.03	23.6 ± 1.73	14.8 ± 0.72	3.16 ± 0.08	21.3 ± 1.31	12.7 ± 0.49
	ExtM	1.29 ± 0.02	6.27 ± 0.43	19.6 ± 1.42	2.31 ± 0.07	6.50 ± 0.29	14.5 ± 0.61	3.98 ± 0.16	7.36 ± 0.31	12.3 ± 0.52
	K-Hype	2.61± 0.07	3.30 ± 0.13	5.40 ± 0.28	3.37 ± 0.11	3.34 ± 0.13	5.13 ± 0.24	3.65 ± 0.13	3.56 ± 0.14	5.11 ± 0.24
	Proposed	1.92 ± 0.04	3.66 ± 0.14	3.21 ± 0.12	3.18 ± 0.11	3.65 ± 0.14	4.99 ± 0.24	3.21 ± 0.11	3.70 ± 0.14	4.95 ± 0.26


Paul Honeine
email: paul.honeine@utt.fr

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

CONSTRAINED REFLECT-THEN-COMBINE METHODS FOR UNMIXING HYPERSPECTRAL DATA

Keywords: Constrained optimization, hyperspectral data, unmixing problem, parallel projection, Cimmino's method

This paper deals with the linear unmixing problem in hyperspectral data processing, and in particular the estimation of the fractional abundances under sum-to-one and nonnegativity constraints. For this purpose, we propose to adapt the reflect-then-combine iterative technique, initially derived by Cimmino. Several strategies are studied in order to handle the constraints, and experimental results are analyzed.

INTRODUCTION

In hyperspectral data processing, spectral unmixing is one of the most fundamental and challenging problems. The spectral unmixing problem consists of breaking down a spectrum into a set of pure spectra, a.k.a. endmembers, and their fractional abundances. Assuming that the endmembers were extracted using any off-the-shelf technique (see [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF] for a survey), the estimation of the abundances brings new opportunities and challenges to both linear [START_REF] Honeine | Geometric unmixing of large hyperspectral images: a barycentric coordinate approach[END_REF] as well as non-linear unmixing problems [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF]. To adopt a physical interpretation, the estimation problem requires the fulfillment of two constraints: the fractional abundances must add up to 100% (a.k.a. sum-toone constraint), and must be additive (a.k.a. non-negativity constraint). See for instance [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF]. This paper deals with constrained fractional abundances estimation, by suitably adapting parallel orthogonal projections/reflections. Methods of orthogonal projections onto subspaces have been bridging the words of geometry and algebra, successfully applied in solving optimization problems. While the alternating projection algorithm goes back at least to Hermann A. Schwarz in the 1870's, it was revisited in the late 1930's by Kaczmarz's cyclic projections [START_REF] Kaczmarz | Angenäherte Auflösung von Systemen linearer Gleichungen[END_REF] and Cimmino's parallel reflections [START_REF] Cimmino | Calcolo approssimato per le soluzioni dei sistemi di equazioni lineari[END_REF]. More recently, there has been a modern revival This work was supported by the French ANR, grant HYPANEMA: ANR-12-BS03-0003.

of projection algorithms, namely with projections onto convex sets within adaptive filtering and machine learning [START_REF] Yukawa | Adaptive filtering based on projection method[END_REF][START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF]. These methods open tremendous opportunities in various domains in signal and image processing.

Hyperspectral data processing still do not take advantage of this increasing research activity, even in the fundamental linear spectral unmixing problem. In this paper, we propose to tackle the major obstacle, which is the fulfillment of the physical interpretation, namely by enforcing non-negativity and sum-to-one constraints within an iterative projection scheme. To this end, we revisit the parallel reflect then combine update rule, initially proposed by Cimmino [START_REF] Cimmino | Calcolo approssimato per le soluzioni dei sistemi di equazioni lineari[END_REF]. We consider two strategies to handle the sum-to-one constraint, either by normalization or by an asymptotic convergence. The nonnegativity constraint is enforced by two different strategies, either by relaxing the reflection or by projecting onto the nonnegative orthant.

THE LINEAR (UN-)MIXING MODEL

Given a spectrum with L wavelength

x = [x 1 x 2 • • • x L ] ⊤ (e.
g. a pixel in a hyperspectral image), the linear mixing model takes the form

x = K k=1 α k m k + ǫ,
where m k is the spectral signature of a pure material, i.e., endmember, and ǫ is the vector of fitness error. In matrix form, we get

x = M α + ǫ, (1)
where

M = [m 1 m 2 • • • m K ], and α = [α 1 α 2 • • • α K ] ⊤
is the vector of abundances to be determined. In this work, it is assumed that the endmembers have been identified, using any off-the-shelf endmember extraction technique. See for instance [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF][START_REF] Honeine | Geometric unmixing of large hyperspectral images: a barycentric coordinate approach[END_REF] and references therein. The (unconstrained) optimization problem, i.e., arg min α x -M α 2 , leads to the optimal solution α = M ⊤ (M M ⊤ ) -1 x, with optimality in the least-squares sense.

In order to adopt a physical interpretation in the unmixing model, the above linear combination must be convex, namely by satisfying the sum-to-one and the non-negativity constraints, respectively 1 ⊤ α = 1 and α ≥ 0, where 1 denotes the unit column-vector of K entries, and the inequality is taken component-wise. An iterative scheme is required to solve this constrained optimization problem, with essentially two steps at each iteration: an approximation step that determines an intermediate solution by minimizing the fitness error; and a step that constrains the solution. The latter step sets to zero the negative values in order to impose the nonnegativity constraint, and normalizes to one to enforce the sum-to-one constraint.

Due to its nature, the sum-to-one constraint can be easily handled in the least-squares optimization problem. To this end, the following augmented model is substituted for (1):

1 x = 1 ⊤ M α + ǫ ǫ . (2)
This model has been recently investigated in the literature for unmixing hyperspectral data, as studied for instance in [START_REF] Honeine | Geometric unmixing of large hyperspectral images: a barycentric coordinate approach[END_REF] using a simplex-based approach. The solution is given by analogy to the above unconstrained problem. Still, the solution is optimal in the least-squares sense, thus leading to a fitness error ǫ within the sum-to-one, as opposed to the explicit constraint enforced in the above iterative scheme. The distinction between these two strategies (explicitly imposed vs. accepting fitness error) is explored in this work in order to derive several constrained reflect-then-combine methods.

CONSTRAINED REFLECT-THEN-COMBINE METHODS

Unconstrained reflect-then-combine method

Back to the unconstrained problem defined in (1), we consider the set of L equations:

x ℓ = m ⊤ ℓ α, for ℓ = 1, 2, . . . , L (3)
where m ⊤ ℓ denotes the vector of the K endmember spectral signatures at the ℓ-th wavelength band. In other words, the ℓ-th row of M is denoted by m ⊤ ℓ . The solution of the unmixing problem is the unique 1 intersection of the L (affine) hyperplanes defined by the above equations.

Cimmino's method consists of two steps to refine the approximate solution α (t) at iteration t, as described by the fol- 1 The intersection is not unique, due to the presence of noise as given by ǫ in (1). However, nothing prevents us from applying the reflection principle, since the affine hyperplane defined by x ℓ = m ⊤ ℓ α + ǫ ℓ is parallel to the one given in the noiseless expression [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF]. lowing update equations:

α (t) /ℓ = α (t) + 2 x ℓ -m ⊤ ℓ α (t) m ℓ 2 m ℓ (4) α (t+1) = L ℓ=1 γ (t) ℓ α (t) /ℓ (5)
The first step involves intermediate solutions, by taking the reflection (i.e., mirror image) of α (t) with respect to each of the hyperplanes defined by [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF]. The second step combines the resulting reflections through weights γ

(t)
ℓ that satisfy the the convexity constraints, i.e., L ℓ=1 γ

(t) ℓ = 1 and γ (t) 1 , γ (t) 2 , . . . , γ (t)
L ≥ 0 at each iteration t. The motivation of this "reflect-then-combine" rule is illustrated in Figure 1. In fact, any convex combination of the reflections provides a better solution that the initial one. In particular, as often applied by Cimmino, one may take the center of gravity of the intermediate solutions, namely by setting γ (t) ℓ = 1/L for all ℓ = 1, 2, . . . , L. It is worth noting that this strategy converges faster than a "project-then-combine" strategy, as shown in Figure 1 (dotted red line vs blue line).

By merging the reflection and combination steps into a single expression, one can write the following update rule:

α (t+1) = α (t) + 2 M ⊤ Γ (x -M α (t) ), (6)
where Γ is a L-by-L diagonal matrix with entries

γ 1 / m 1 2 , γ 2 / m 2 2 , . . . , γ L / m L 2 .
This method has several desirable properties, such as parallelization as given by the update rule (4) which can be processed separately on a L-core machine, as well as multiple pixel processing. For the latter case, by collecting the x and the corresponding α (t) of each pixel into X and A (t) respectively, we get the following update rule:

A (t+1) = A (t) + 2 M ⊤ Γ (X -M A (t) ).

Sum-to-one constraint

Next, we explore two strategies to handle the sum-to-one constraint, namely K k=1 α (t+1) k = 1 at each iteration t. In the first strategy, the system converges to this constraint by using the augmented model, while in the second one its is verified by normalization each iteration.

In the first strategy, we consider the augmented model (2), which can be easily investigated using the "reflect-thencombine" update rule. The resulting rule can be written as

α (t+1) = α (t) + 2 1 M ⊤ Γ 1 x - 1 ⊤ M α (t) ,
[S augment ] where Γ is here a (L + 1)-by-(L + 1) diagonal matrix. In this paper, the augmented model extends naturally the "reflectthen-combine" update rule. This strategy can be viewed as

x ℓ = m ⊤ ℓ α xℓ ′ = m ⊤ ℓ ′ α α (t) α (t) /ℓ α (t) /ℓ ′ α (t+1)
←optimal solution Fig. 1. Illustration in two-dimensions of the "reflect-then-combine" update rule, with two hyperplanes (black lines). The final solution α (t+1) (⋆) is obtained by a convex combination of the reflections (•) of the initial solution α (t) (⋆). As shown, any combination of the reflections (-) is a better solution that every combination of the projections (• • • ).

including a new reflection at each iteration, with respect to the (affine) hyperplane defined by the equation 1 = 1 ⊤ α (t) .

In the second strategy, we propose to operate a normalization after every adaptation, namely by applying successively the following rules at each iteration t:

α ( * ) = α (t) + 2 M ⊤ Γ (x -M α (t) ), α (t+1) = α ( * ) 1 ⊤ α ( * ) , [S normalize ]
where α ( * ) is the intermediate (unnormalized) solution. One can also view the above normalization step as a projection, onto the (affine) hyperplane defined by the equation 1 = 1 ⊤ α (t) . But this is not an orthogonal projection, and therefore it does no longer satisfy the non-expansivity property. This means that the distances between two projected entries are not guarantied to be less than the distances of the original ones. While this may affect the convergence, the combination of the reflection and projection seems fortunately to work well.

There exists two major differences between these two strategies: On the one hand, the normalization is a nonorthogonal projection onto the hyperplane, as opposed to the orthogonal reflection with the augmented model. On the other hand, normalization guarantees the sum-to-one constraint at each iteration, as opposed to an asymptotic convergence with the previous approach.

Non-negativity constraint

Next, we derive two strategies to enforce the non-negativity constraint, namely α ≥ 0 at each it-eration t. The first one relaxes the reflections to prevent getting out of the non-negative orthant, while the second one operates projection onto it.

In the first strategy, we revisit the reflection step given in the update rule (4), by including a relaxation with

α (t) /ℓ = α (t) + 2 η ℓ x ℓ -m ⊤ ℓ α (t) m ℓ 2 m ℓ , [S relax ]
where the relaxation weight η ℓ ∈ [0 ; 1] is chosen in order to impose the non-negativity of entries in α

/ℓ . In practice, it is either set to 1 when the reflection satisfies the constraint; otherwise it is chosen on the boundary, namely when one of the entries of α (t) /ℓ goes to zero. In other words,

η ℓ = min k -m ℓ 2 x ℓ -m ⊤ ℓ α (t) [α (t) ] k [ m ℓ ] k when the right-hand-side is within [0 ; 1]; otherwise η ℓ = 1.
The second strategy consists of setting to zero all negative values after each reflect-and-combine. This corresponds to apply the following rules at each iteration t:

α ( * ) =α (t) + 2 M ⊤ Γ (x -M α (t) ), α (t+1) =[α (t+1) 1 • • • α (t+1) K ] ⊤ with α (t+1) k = max{α ( * ) k ; 0} [S set-to-0 ]
where α ( * ) is the intermediate (unconstrained) solution.

A major difference between these two strategies is that the latter enforces null values to prevent negative ones. The former strategy leads to non-null values, due to the combination process.

Fully-constrained optimization

In the previous sections, we have proposed 2 strategies to handle the sum-to-one constraint, and 2 strategies to enforce the non-negativity constraint. This leads to 4 variants for the fully-constrained optimization problem. The following table summarizes the properties of these fully-constrained algorithms:

sum-to-one non-negativity

S augment & relax 1 ⊤ α → 1 α > 0 S augment & set-to-0 1 ⊤ α → 1 α ≥ 0 S normalize & relax 1 ⊤ α = 1 α > 0 S normalize & set-to-0 1 ⊤ α = 1 α ≥ 0
This table illustrates that the hard constraints of sum-toone and non-negativity are replaced by soft constraints, i.e., convergence. The price to pay for fulfilling the hard constraints, as given by the S normalize & set-to-0 strategy, is a slower slower convergence rate. This is essentially due to the normalization in S normalize , which corresponds to a non-orthogonal projection and therefore loosing the non-expansivity property. Fig. 2. Convergence of the fully-constrained optimization methods: S augment & relax ("-"), S augment & set-to-0 ("-") and S normalize & relax ("--") strategies. The dotted lines correspond to the real abundances α = [0.4 0.6 -0.05] ⊤ . The estimates using non-negative least-squares (+) and the fullyconstrained least-squares ( * ) are also given.

In this paper, we consider also the soft constraints. It is worth noting that most of the work on hyperspectral unmixing problem consider relaxing the constraints. See for instance [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF].

EXPERIMENTATIONS

To avoid the dilemma of the lack of ground truth information for comparing unmixing methods, and due to space limitation, we considered the linear combination of real hyperspectral signatures extracted from the USGC library: grass, cedar and asphalt. These spectra consist of 2151 bands covering wavelengths ranging from 0.35 to 2.5 µm. A spectral mixture was generated by a linear combination of three pure spectral signatures with fractional abundances α = [0.4 0.6 -0.05] ⊤ . The data were corrupted by an additive Gaussian noise with a SNR ≈ 35 dB.

The convergence of the fully-constrained optimization methods is illustrated in Figure 2 for 100 iterations. As expected, the augmented model as in S augment is preferred to the normalization S normalize . Both the S augment & relax strategy ("-") and the S augment & set-to-0 strategy ("-") converge faster than the S normalize & relax strategy ("--"). The S normalize & set-to-0 strategy was discarded due to its poor (and bad) convergence properties. For a comparative study, we also included in the figure the estimates obtained from the non-negative least-squares method ("+") and the fully-constrained least-squares method (" * "). See for instance [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF].

CONCLUSION AND PERSPECTIVES

This work shows that conventional Cimmino's reflect-thencombine iterative scheme can be adapted for solving the constrained unmixing problem. Several strategies are studied in order to handle the sum-to-one and non-negativity constraints, and experimental results are analyzed.

As for future work, we are studying the optimal choice of the weights in the combination stage. We seek also to optimize its selection in order to satisfy both sum-to-one and non-negativity constraints. In the latter case for instance, one may allow negative values in reflections as long as there exists a non-negative convex combination. We are also conducting an analysis of the algorithm, in the same spirit as in [START_REF] Chen | Nonnegative least-mean-square algorithm[END_REF]. It is also interesting to study Kaczmarz's cyclic projections, as well as random projections/reflections. It is clear that this first work opens the way to a new class of algorithms in adaptive filtering and machine learning [START_REF] Yukawa | Adaptive filtering based on projection method[END_REF][START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF].


Paul Honeine
email: paul.honeine@utt.fr

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

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

CONSTRAINED KACZMARZ'S CYCLIC PROJECTIONS FOR UNMIXING HYPERSPECTRAL DATA

Keywords: Constrained optimization, Kaczmarz's cyclic projections, hyperspectral data, unmixing problem

The estimation of fractional abundances under physical constraints is a fundamental problem in hyperspectral data processing. In this paper, we propose to adapt Kaczmarz's cyclic projections to solve this problem. The main contribution of this work is two-fold: On the one hand, we show that the non-negativity and the sum-to-one constraints can be easily imposed in Kaczmarz's cyclic projections, and on the second hand, we illustrate that these constraints are advantageous in the convergence behavior of the algorithm. To this end, we derive theoretical results on the convergence performance, both in the noiseless case and in the case of noisy data. Experimental results show the relevance of the proposed method.

INTRODUCTION

Due to low spatial resolution of hyperspectral cameras (i.e., imaging spectrometers), acquired spectra are mixtures of spectra of some pure materials. These pure components are also called endmembers. The unmixing problem consists of breaking down a spectrum into the pure spectra of endmembers and their fractional abundances. A wide range of geometrical and statistical methods has been devised for the unmixing problem. See [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF] for an extensive overview. Hyperspectral data processing often suffers from a number of limitations, which includes computation complexity and the ease of implementation, making iterative and parallel algorithms suitable for implementation [START_REF] Heylen | Fully constrained least-squares spectral unmixing by simplex projection[END_REF][START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF].

Assuming that the endmembers were extracted using any off-the-shelf technique, the problem of estimating the fractional abundances provides new opportunities and challenges to both linear [START_REF] Honeine | Geometric unmixing of large hyperspectral images: a barycentric coordinate approach[END_REF] as well as non-linear unmixing problems [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF]. In order to have a coherent physical interpretation, two hard constraints need to be satisfied in the estimation problem: the This work was supported by the French ANR, grant HYPANEMA: ANR-12-BS03-0003.

sum-to-one of all the fractional abundances and there nonnegativities. In this work, we study the problem of estimating constrained fractional abundances, by suitably adapting a cyclic orthogonal projection scheme.

Methods of orthogonal projections onto spaces, and more recently onto convex sets, have been successfully applied for solving many optimization problems, in the field of signal processing in particular. These methods include the celebrated Cimmino's parallel projections [START_REF] Cimmino | Calcolo approssimato per le soluzioni dei sistemi di equazioni lineari[END_REF] and Kaczmarz's cyclic projections [START_REF] Kaczmarz | Angenäherte Auflösung von Systemen linearer Gleichungen[END_REF]. More recently, projections onto spaces have been generalized to projections onto convex sets. These methods provide a unified framework for tackling many signal processing problems, including adaptive filtering and machine learning [START_REF] Yukawa | Adaptive filtering based on projection method[END_REF][START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF]. To the best of our knowledge, these methods have not yet been exploited in hyperspectral image processing. Very recently, we have successfully adapted Cimmino's parallel projections [START_REF] Cimmino | Calcolo approssimato per le soluzioni dei sistemi di equazioni lineari[END_REF] for the hyperspectral unmixing problem [START_REF] Honeine | Constrained reflect-then-combine methods for unmixing hyperspectral data[END_REF].

Coupled with these advances, this paper shows that the spectral unmixing problem can take advantage of these developments. To this end, we propose to adapt Kaczmarz's cyclic projection for estimating the fractional abundances. Appropriate strategies for enforcing the non-negativity and sum-toone constraints are advanced as a natural extension of the conventional unconstrained Kaczmarz's method: the sum-to-one constraint is enforced by a projection scheme, while the nonnegativity is satisfied by relaxing the projections.

THE LINEAR UNMIXING MODEL

Given a spectrum with L wavelength

x = [x 1 x 2 • • • x L ]
(e.g. a pixel in a hyperspectral image), the linear mixing model takes the form

Mα = x + , ( 1
)
where M is the L-by-K matrix of the K spectral signatures of a pure material, i.e., endmember, and is the vector of fitness error, and α is the vector of K abundances to be determined. In this work, we assume that the endmembers have been properly identified, using a supervised or an unsupervised endmember extraction technique. The paper [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF] provides an extensive overview of existing techniques and stateof-the-art methods.

To adopt a physical interpretation of the unmixing problem, the linear combination in model (1) needs to be convex. In other words, the following constraints are required:

• The sum-to-one of the contributions, 1 α = 1, where 1 denotes the unit column-vector of K entries;

• and the non-negativity of the contributions, α ≥ 0, where the inequality is taken component-wise.

An iterative scheme is required to solve the constrained optimization problem: arg min α Mαx 2 , under the above constraints. See for instance [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF].

CONSTRAINED KACZMARZ'S CYCLIC PROJECTIONS

To revisit the unconstrained Kaczmarz's method, we consider the noiseless model, defined by

M α = x. (2)
In other words, Kaczmarz's method assumes that x belongs to the range of the matrix M . The impact of the presence of noise, as given in model [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF], is studied in next Section.

Unconstrained Kaczmarz's cyclic projections

The linear system (2) of L equations and K unknowns gives rise to a set of L affine hyperplanes,

H = α m α = x . (3)
In this expression, m denotes the -th row of M , namely the spectral signatures of the endmembers at the -th wavelength band. The solution of the linear system (2) is the intersection of the L affine hyperplanes, i.e., H 1 ∩H 2 ∩. . .∩H L , where the intersection is a unique point in the noiseless case. It is worth noting that M is full column rank by construction, since non of the endmembers can be written as a linear combination of other endmembers. The conventional Kaczmarz's method sweeps through the affine hyperplanes, by projecting orthogonally onto each and taking this as the next iterate. In a mathematical form, at the tth iteration, the projection of α t-1 onto the affine hyperplane H reads

α t-1⊥ = α t-1 + x -m α t-1 m 2 m , (4)
with α 0 some initial guess. This process is iterated until a given convergence criterion is satisfied. In this expression, is selected at each iteration t, i.e., a slight abuse of notation is made in using rather than (t).

1 H m H α t-1 α t-1⊥ α t θ α opt
Several strategies [START_REF] Needell | Randomized Kaczmarz solver for noisy linear systems[END_REF][START_REF] Strohmer | A randomized Kaczmarz algorithm with exponential convergence[END_REF][START_REF] Censor | A note on the behavior of the randomized Kaczmarz algorithm of Strohmer and Vershynin[END_REF] have been proposed to choose the sequence of selected affine hyperplanes, i.e., the choice at iteration t of the -th row in the matrix M . The classical cyclic selection, i.e., = t mod L, is the simplest one; however it could lead to low convergence. Other strategies have been proposed to overcome this drawback, the most known are the maximum error selection and the randomization. The maximum error selection requires the evaluation of the numerator in ( 4) for all at each iteration t, which leads to a high computational burden. In randomization is chosen randomly with probability proportional to the denominator in (4), i.e., m 2 . This strategy has shown very interesting properties [START_REF] Strohmer | A randomized Kaczmarz algorithm with exponential convergence[END_REF], but is also criticized in [START_REF] Censor | A note on the behavior of the randomized Kaczmarz algorithm of Strohmer and Vershynin[END_REF].

This paper shows that, in the case of a constrained solution as in the unmixing problem, Kaczmarz's method does not suffer from the problem of appropriate selection of the affine hyperplanes. Moreover, we provide theoretical results on the convergence, for the noiseless and the noisy cases.

Sum-to-one constraint

To impose the sum-to-one constraint, we consider a two-step strategy at each iteration t: a projection of α t-1 onto the affine hyperspace H , yielding α t-1⊥ , followed by a normalization to satisfy the constraint, yielding α t . The proposed strategy is illustrated in Figure 1.

Let H be the affine hyperplane defined by the sum-to-one, namely H = α | 1 α = 1 . By analogy with the expression of the projection (4), the normalization of any α t-1⊥ is given by

α t = α t-1⊥ + 1 K 1 -1 α t-1⊥ 1, (5)
where m is now substituted with 1, the unit column-vector of K ones. By combining this projection [START_REF] Cimmino | Calcolo approssimato per le soluzioni dei sistemi di equazioni lineari[END_REF] with the projection onto H , as given in (4), at iteration t, we get:

α t = α t-1⊥ + 1 K 1 -1 α t-1+⊥ 1 = α t-1 + x -m α t-1 m 2 m + 1 K 1 -1 α t-1 - x -m α t-1 m 2 1 m 1.
Since we have the identity 1 α t-1 = 1 from the previous iteration, we obtain the following update rule:

α t = α t-1 + x -m α t-1 m 2 I -1 K 11 m . ( 6
)
Here, I denotes the matrix identity of size K-by-K, and the matrix I -1 K 11 is often called centering matrix in the machine learning literature.

It is worth noting that there is another way to satisfy the sum-to-one constraint, by normalizing using the form:

α t = α t-1⊥ 1 α t-1⊥ .
This normalization can also be viewed as a non-orthogonal projection, and therefore looses the property of nonexpansive transformation. Technically, this means that the distance between two projections may exceed the distance between the two original elements, and therefore compromising the speed of convergence. For these reasons, this normalization is not considered in this paper, as we adopt the (orthogonal) projection given in (5).

Imposing also the non-negativity constraint

To enforce the non-negativity constraint, namely α t ≥ 0 at each iteration, we propose to relax the projection in [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF], by considering the following update rule:

α t = α t-1 + η t x -m α t-1 m 2 I -1 K 11 m . (7)
The relaxation parameter η t ∈ [0 ; 1] is chosen in order to impose the non-negativity of entries in α t .

To this end, let

η t = min i=1,...,K η (i) t ,
where η

(i)
t is a valid stepsize in the i-th direction. The value of the latter depends on one of the two cases:

• If x -m αt-1 m 2 I -1 K 11 m i > 0,
there is no constraint on the stepsize, and therefore could be set to η

(i) t = 1;
• Otherwise, the stepsize should be reduced such that

η (i) t ≤ [αt-1]i I- 1 K 11 m i .
In order to account for noise, we consider always a relaxation of the projection. The stepsize is set to a smaller value than one (e.g., 0.1 or below depending on the above relations), as recommended in adaptive filtering literature. Essentially, the update rule ( 7) is a normalized least mean squares (NLMS) algorithm applied on centered data. Other strategies to enforce the non-negativity constraint can also be applied, such as a post-processing by replacing every negative entry in α t with zero. However, this strategy of projection onto the positive half-space may destroy the sumto-one constraint, and therefore it is not recommended.

THEORETICAL RESULTS

In this section, we provide theoretical results, namely on the convergence rate in both the noiseless (2) and the noisy cases (1). To this end, we consider the influence of both sum-to-one and non-negativity constraints by using the update rule [START_REF] Heylen | Fully constrained least-squares spectral unmixing by simplex projection[END_REF].

First of all, the impact of the non-negativity constraint is clear. All the results are enforced to be on the non-negative orthant in the K-dimensional space. Furthermore, the studied hyperspectral unmixing problem has the following interesting property: Since m is a vector of values from spectral signatures at the -th wavelength band, we have m ≥ 0, and therefore all affine hyperplanes defined by (3) have nonnegative normal vectors. Moreover, the sum-to-one affine hyperplane H is the bisector of the orthant, and therefore the angle between it and any H cannot exceed π/4. Let θ be the angle between the two affine hyperplanes H and H, then cos(θ ) > 1/ √ 2, for all = 1, 2, . . . , L. Let α opt be the (unknown) optimal solution of the noiseless problem (2).

The noiseless case

The following theorem provides the rate of convergence for the proposed constrained Kaczmarz' method, in the case of the noiseless model [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]. To this end, we consider the influence of the sum-to-one constraint on the convergence.

Theorem 1 (Convergence in the noiseless case).

Starting from an initial guess α 0 , the algorithm converges to the optimal solution α opt at the rate

α t -α opt = α 0 -α opt (t ) t ≤t cos 2 (θ ),
where θ is the angle between H and H, and the product is taken over all visited affine hyperplanes until iteration t.

Proof. As in the conventional Kaczmarz's method, it is assumed that the solution of the linear system exists, is unique, and belongs to the set H . On the one hand, we have from projection (i.e., according to the Pythagorean theorem)

cos(θ ) = α t-1⊥ -α opt α t-1 -α opt ,
where cos(θ ) is positive as discussed above. On the other hand, we also have

cos(θ ) = α t -α opt α t-1⊥ -α opt .
By combining both expressions, we get

α t -α opt = α t-1⊥ -α opt cos(θ ) = α t-1 -α opt cos 2 (θ ) . . . = α 0 -α opt (t ) t ≤t cos 2 (θ ).
This concludes the proof. This theorem demonstrates the rate of convergence of the proposed algorithm. It is worth noting that cos(θ ) can be easily evaluated, since cos(θ ) = 1 m 1 m , where the numerator is the sum of the entries of m , and the denominator involves the sum of their squares.

Moreover, each iteration in the proposed algorithm is essentially a double projection. This alternating projection can be applied infinitely on a fixed , i.e., affine hyperplane H , without the need to sweep through all the affine hyperplanes. In this case, the Theorem 1 becomes for a fixed : α tα opt = α 0α opt (cos(θ ))

2t . In practice, due to the presence of noise, it is not enough to use a single , but rather to sweep over all the hyperplanes. The next Section provides a study of the impact of noise on the proposed algorithm.

The noisy case

Very few works consider the influence of noise in the conventional (unconstrained) Kaczmarz's cyclic projections. The presence of noise is seldom considered; see [START_REF] Needell | Randomized Kaczmarz solver for noisy linear systems[END_REF] for the randomized Kaczmarz's method. In this section, we provide a theoretical analysis of the presence of an error in the proposed constrained optimization problem. We assume that the endmember spectra have been well estimated, and the noise is on the investigated spectrum as given in (1), namely Mα = x + , where = [ 1 2 . . . L ] is the error vector. The solution of this noisy problem belongs to the intersection of the L affine hyperplanes

H 1 ∩ H 2 ∩ . . . ∩ H L , defined by H = α m α = x + . H H H α t-1 α t-1⊥ α t-1⊥ α t α t m m 2 θ θ α opt
Fig. 2. Illustration in two-dimensions of the presence of noise. In this case, H is not available (neither all gray elements), but the affine hyperplane H associated to the noisy data. Therefore, at instant t-1, the approximation α t-1 leads to α t , rather than α t .

It is easy to see that these affine hyperplanes can also be defined from the hyperplanes associated for the noiseless case, with (See for instance [START_REF] Needell | Randomized Kaczmarz solver for noisy linear systems[END_REF]Lemma 2.2]):

H = α + m 2 m α ∈ H . (8)
Let α t-1 be the estimate obtained at iteration t in presence of noise, and α t-1⊥ its projection onto H . The following theorem is a generalization of Theorem 1 to the case of noisy data. See Figure 2 for an illustration.

Theorem 2 (Convergence in the noisy case). Starting from an initial guess α 0 , the algorithm converges at the rate

α t -α opt 2 = α 0 -α opt 2 (t ) t ≤t cos 4 (θ ) + (t ) t ≤t 2 m 2 cos 2 (θ ) (t ) t <t cos 4 (θ ),
where θ is the angle between H and H (as well as H ).

Proof. First, we have from projection the following relation: cos(θ ) = α tα opt / α t-1⊥α opt . By applying the Pythagorean theorem, the above denominator can be decomposed into

α t-1⊥ -α opt 2 = α t-1⊥ -α opt 2 + m 2 m 2 ,
where the orthogonality and the definition (8) are used. The second term in the right-hand-side is simply 2 / m 2 . The first term in the right-hand-side is (see the proof of Theorem 1): α t-1⊥α opt = α t-1α opt cos(θ ). By combining these expressions, we get

α t -α opt 2 = α t-1 -α opt 2 cos 4 (θ ) + 2 m 2 cos 2 (θ ) . . . = α 0 -α opt 2 (t ) t ≤t cos 4 (θ ) + (t ) t ≤t 2 m 2 cos 2 (θ ) (t ) t <t cos 4 (θ ).
This concludes the proof. This theorem demonstrates the impact of the selection, i.e., (t) at each iteration t, on the approximation error. In fact, thanks to the geometric progression of the cosines, the hyperplanes considered in the first iterations have their errors highly reduced, as opposed to the last projections with less "cosine" weights. In other words, it is clear that the bandwidths in the hyperspectral data with the largest errors, i.e., large values of 2 / m 2 , should be used in the first iterations. Therefore, this theorem provides a new selection criterion for the cyclic projections, optimal for the studied constrained optimization problem.

EXPERIMENTATIONS

The proposed learning algorithm was used for hyperspectral data unmixing, generated by a linear combination of three pure materials with abundance α = [0.4 0.65 -0.05] , where we injected a sign error on the third abundance. These materials are grass, cedar and asphalt, with spectral signatures extracted from the USGC library. These spectra consist of 2151 bands covering wavelengths ranging from 0.35 to 2.5 µm. The data were corrupted by an additive Gaussian noise. Figure 3 illustrates the convergence of the proposed algorithm, for a SNR of ≈ 25 dB and ≈ 18 dB, with the largest value of the stepsize fixed to 0.1 and a single sweep over all the bandwidths. For a comparative results, we also show in the same figure the estimates using non-negative least-squares and the fully-constrained least-squares techniques [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF].

CONCLUSION AND PERSPECTIVES

In this paper, we presented the development of a new class of unmixing methods, by providing constrained counterpart of projection-based algorithms [START_REF] Yukawa | Adaptive filtering based on projection method[END_REF][START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF]. New theoretical results were derived for the studied constrained optimization problem, including the impact of noise on the analyzed spectrum. Experimentations collaborate these results. As for future work, we seek to study the impact of noise on the endmembers. We are also conducting a statistical analysis on the impact of noise, as well as connections with adaptive filtering. The estimates using non-negative least-squares (given by "+") and the fully-constrained least-squares (given by " ") are also given.

Fig. 1 .

 1 Fig.1. Illustration in two-dimensions of the constrained Kaczmarz's cyclic projections. The affine hyperplane of the sum-to-one constraint is H (black line), i.e., the hyperplane defined by the normal vector 1.

Fig. 3 .

 3 Fig. 3. Convergence of the constrained Kaczmarz's optimization method, for a mixture of two spectra with additive noise. The dotted lines correspond to the real abundances α = [0.4 0.65 -0.05] .The estimates using non-negative least-squares (given by "+") and the fully-constrained least-squares (given by " ") are also given.