Zineb Noumir

Paul Honeine

Cédric Richard

ONE-CLASS MACHINES BASED ON THE COHERENCE CRITERION

Keywords: support vector machines, machine learning, kernel methods, one-class classification

The one-class classification problem is often addressed by solving a constrained quadratic optimization problem, in the same spirit as support vector machines. In this paper, we derive a novel one-class classification approach, by investigating an original sparsification criterion. This criterion, known as the coherence criterion, is based on a fundamental quantity that describes the behavior of dictionaries in sparse approximation problems. The proposed framework allows us to derive new theoretical results. We associate the coherence criterion with a one-class classification algorithm by solving a least-squares optimization problem. We also provide an adaptive updating scheme. Experiments are conducted on real datasets and time series, illustrating the relevance of our approach to existing methods in both accuracy and computational efficiency.

INTRODUCTION

In machine learning, several problems exhibit only a unique class for training. The decision rule consists of identifying if a new instance belongs to the learnt class or to an unknown class. This machine learning problem is the so-called oneclass classification problem [START_REF] Tax | One-class classification[END_REF][START_REF] Schölkopf | Estimating the support of a high-dimensional distribution[END_REF]. It has been applied with success for novelty detection, and extends naturally to tackle multiclass tasks by learning each class separately [START_REF] Zhang | Adaptive and online one-class support vector machine-based outlier detection techniques for wireless sensor networks[END_REF][START_REF] Ratle | A comparison of one-class classifiers for novelty detection in forensic case data[END_REF].

Several methods have been developed to solve this problem, the most widely studied being the one-class support vector machines (SVM). It determines a sphere of minimum volume that encloses all (or most of) the training data, by estimating its center and radius. One-class SVM takes advantage of many properties from SVM literature, such as the nonlinear extension by using kernel functions and the sparse solution. The sparsity property states that the center of the sphere explores only a small fraction of the training samples, known as support vectors (SVs). A quadratic programming technique is often applied to solve this problem, i.e., identifying the SVs and estimating the optimal parameters. This work was supported by ANR-08-SECU-013-02 VigiRes'Eau.

In this paper, we derive a one-class classification machine based on a new sparsification rule, the coherence criterion. This criterion is based on the coherence parameter, a fundamental quantity that describes the behavior of dictionaries in sparse approximation problems. We have recently investigated this sparsification criterion, and applied it with success in nonlinear filtering and online prediction in time series. See [START_REF] Richard | Online prediction of time series data with kernels[END_REF][START_REF] Honeine | On-line nonlinear sparse approximation of functions[END_REF]. The coherence criterion provides an elegant model reduction criterion with a low computational demanding procedure. This framework allows us to derive new theoretical results, mainly an upper bound on the error of approximating the center by the resulting reduced model. We associate the coherence criterion with a one-class classification algorithm by solving a least-squares optimization problem. We also provide an adaptive updating scheme for incrementing or decrementing the model order.

The rest of the paper is organized as follows. Section 2 outlines the classical one-class SVM. We describe our approach in Section 3, by providing theoretical results and deriving an appropriate one-class classification algorithm. Section 4 provides an experimental study on real datasets and time series.

ONE-CLASS SVM

Let x 1 , x 2 , . . . , x n be the set of training samples, et let Φ(•) be a nonlinear transformation defined by the use of a kernel, κ(x i , x j ) = Φ(x i ), Φ(x j ) . The one-class SVM finds a sphere, of minimum volume, containing all (or most of) the training samples. This sphere, described by its center c and its radius r, is obtained by solving the constrained optimization problem:

min r,c,ζ r 2 + 1 νn n i=1 ζ i subject to Φ(x i ) -c 2 ≤ r 2 + ζ i for all i = 1, 2, . . . , n
In this expression, ζ 1 , ζ 2 , . . . , ζ n is a set of non-negative slack variables and ν a positive parameter that specifies the tradeoff between the sphere volume and the number of outliers, i.e., samples lying outside the sphere. By introducing the Karush-Kuhn-Tucker optimality conditions, we get

c = n i=1 α i Φ(x i ), (1)
where the α i 's are the solution to the optimization problem:

max α n i=1 α i κ(x i , x i ) - n i,j=1 α i α j κ(x i , x j ) subject to n i=1 α i = 1 and 0 ≤ α i ≤ 1 νn for all i = 1, 2, . . . , n.
It is well known that only a small fraction of the training samples contributes to the model (1). These samples, called support vectors (SVs), correspond to data outside or lying on the boundary of the sphere, and consequently contribute to the definition of the radius r. Thus, any new sample x that satisfies Φ(x)c > r can be considered as an outlier, where the distance is

Φ(x)-c 2 = i,j∈I α i α j κ(x i , x j )-2 i∈I α i κ(x i , x)+κ(x, x),
where I denotes the set of SVs, i.e., training samples with non-zero α i 's. Let |I| denotes its cardinality.

PROPOSED ONE-CLASS METHOD

Let c n denotes the center of the n samples, namely

c n = 1 n n i=1 Φ(x i ). (2)
We consider a sparse model of the center with

c I = i∈I α i Φ(x i ), (3)
where the problem consists of properly identifying the subset I ⊂ {1, 2, . . . , n} and estimating the optimal coefficients α i 's. While classical one-class SVM operates a joint optimization, we consider in this paper a separate optimization scheme:

1. Identify the most relevant samples in the expansion (3), by using the coherence parameter in the sparsification rule.

2. Estimate the optimal coefficients, with optimality in the least-squares sense, namely by minimizing c n -c I 2 .

Sparsification rule with the coherence criterion

The coherence of a set {Φ(x i ) | i ∈ I} is defined by the largest absolute value of the off-diagonal entries of the Gram (kernel) matrix, namely a µ-coherent set is1

µ = max i,j∈I i =j |κ(x i , x j )|.
With the model order being fixed in advance, we consider the set of least coherence as the relevant set in the expansion (3). We have previously applied with success this coherence criterion in nonlinear filtering. See [START_REF] Richard | Online prediction of time series data with kernels[END_REF][START_REF] Honeine | On-line nonlinear sparse approximation of functions[END_REF]. In order to consider this sparsification rule for the one-class problem, one needs to study the relevance of the sparse model c I (obtained by the coherence criterion) with respect to the full-order center c n in (2). The following theorem provides a guarantee that this approximation error is upper bounded.

Theorem 1. For the sparse solution c I satisfying the coherence criterion, the approximation error can be upper bounded with

c n -c I ≤ 1 -|I|/n max i κ(x i , x i ) -µ. ( 4
)
This upper bound takes the form (1 -|I|/n) √ 1 -µ for unitnorm kernels.

Proof. Let P I be the projection operator onto the space spanned by the elements Φ(x i ) for i ∈ I, thus c I = P I c n . Then, we have

c n -c I = 1 n n i=1 (1 -P I )Φ(x i ) ≤ n i=1 1 n Φ(x i ) -P I Φ(x i ) = 1 n i ∈I Φ(x i ) -P I Φ(x i )
where the inequality is due to the generalized triangular inequality, and the last equality to the fact that Φ(x i ) -P I Φ(x i ) = 0 for all Φ(x i ) belonging to the expansion in c I , i.e., for i ∈ I. Moreover, we have

Φ(x i ) -P I Φ(x i ) 2 = Φ(x i ) 2 -PΦ(x i ) 2 = κ(x i , x i ) -max γ j∈I γ j κ(x j , x i ) j∈I γ j Φ(x j ) ≤ κ(x i , x i ) -max k∈I |κ(x k , x i )| κ(x k , x k ) ≤ κ(x i , x i ) -µ
where the first equality is due to the Pythagorian theorem, and the second equality follows the fact that the square norm of the projection of Φ(x i ) corresponds to the maximum

Table 1. Experimental results with the classification error for each one-class classifier, and the mean error, as well as the total computational time (in seconds). The ratio of common SVs, with respect to the results obtained from the classical one-class SVM, for each of the proposed methods is given, as well as the mean ratio of common SVs. scalar product Φ(x i ), ϕ over all the unit-norm functions ϕ, namely ϕ = j∈I γ j Φ(x j )/ j∈I γ j Φ(x j ) . The first inequality results from a specific distribution of the coefficients, with γ j = 0 for all j ∈ I except for a single index k with γ k = ±1, depending on the sign of κ(x k , x i ). The last inequality follows from the coherence criterion.

Optimal parameters

The coefficients are estimated by minimizing the approximation error c nc I , namely

α = arg min αi 1 n n i=1 Φ(x i ) - i∈I α i Φ(x i ) 2 ,
where α is a column vector of the optimal coefficients α k 's for k ∈ I. Taking the derivative of this cost function with respect to each α k , and setting it to zero, we get

1 n n i=1 κ(x k , x i ) = i∈I α i κ(x k , x i ), for every k ∈ I.
In matrix form, we obtain

α = K -1 κ, ( 5
)
where K is the Gram (kernel) matrix, with entries κ(x i , x j ) for i, j ∈ I and κ is a column vector with entries

1 n n i=1 κ(x k , x i ) for k ∈ I.
An interesting property of the proposed method is that one can easily increment or decrement the model order, without the need to operate from scratch as in SVM. For instance2 , one does not need to inversion the new Gram matrix. Let K -1 m be the Gram matrix defined on the set I with cardinality m. If one wishes to increment the model order to m + 1 by injecting the entry x k , then we have

K m+1 = K m b b κ(x k , x k )
where b is the vector with entries κ(x k , x i ) for all i ∈ I.

Then, the inverse of K m+1 can be computed using the block matrix inversion identity, with

K -1 m+1 = K -1 m 0 0 0 + 1 c -K -1 m b 1 -b K -1 m 1 , (6)
where c = κ(x k , x k )b K -1 m b and 0 is a m-length column vector of zeros. The decremental update can be easily derived since, by using the notation

K -1 m =
Q m-1 q q q 0 , we obtain from (6) the update K -1 m-1 = Q m-1qq /q 0 .

EXPERIMENTATION

Multiclass classification

We have tested on three real datasets from the UCI machine learning repository, and well-known in the literature of oneclass machines [START_REF] Wang | Structured one class classification[END_REF]: "iris" with 150 samples in 3 classes and 4 features (only third and fourth features are often investigated), "wine" with 178 samples in 3 classes and 13 features, and "breast cancer" from Wisconsin with 683 samples in 2 classes and 9 features.

Fig. 1 .

 1 Fig. 1. Time series of the chlorine concentration, with one day for training (blue) and the next 3 days for error estimation (red).

This definition corresponds to a unit-norm kernel, i.e., κ(x, x) = 1 for every x; otherwise, substitute κ(x i , x j )/ κ(x i , x i ) κ(x j , x j ) for κ(x i , x j ) in the expression.

Updating expressions for the vectors α and κ can be easily derived, and were omitted due to space limitation.

In order to provide a comparative study, we considered the same configuration as in [START_REF] Liu | Fast support vector data descriptions for novelty detection[END_REF]. The Gaussian kernel was applied, with κ(x i , x j ) = exp(x ix j 2 /2σ 2 ) where σ is the tunable bandwidth parameter. For an -class classification task, a one-class classifier was constructed for each class, called target class. The target set was randomly partitioned into two subsets, one used for training and the other for the target test set. The whole test set contains the target test set and all the samples from the other classes. The classification error were estimated with the optimal parameters obtained by a ten-fold cross-validation on a grid search over ν ∈ {2 -5 ; 2 -4 ; • • • ; 2 4 ; 2 -1.5 } and σ ∈ {2 -5 ; 2 -4 ; • • • ; 2 4 ; 2 5 }.

To compare the proposed method to one-class SVM (which is also comparable to [START_REF] Liu | Fast support vector data descriptions for novelty detection[END_REF]), we first considered the same number of SVs as in SVM. Table 1 gives the classification errors for both methods, as well as the ratio of shared SVs. We found that our approach is very competitive with the same model order as SVM, with up to 80% of common SVs. We also conducted a series of experiments to identify the optimal number of SVs for our algorithm, selected from {3, 4, . . . , 15}. This yields significant accuracy, as given in Table 1. The computational cost, using the best configuration of the tunable parameters, are given in terms of CPU time 3 .

Time series domain description

We also conducted some experiments on a time series. It consists of the variation in chlorine concentration at a given node in a water network. Chlorine is a highly efficient disinfectant, injected in water supplies to kill residual bacteria. This time series exhibits large fluctuations due to the variations in water consumption and an inefficient control system. This data, taken from the public water supplies of the Cannes city in France, was sampled at the rate of a sample every 3 minutes. We considered 4 days of chlorine concentration measures. See Figure 1.

To capture the structure of the time series, a 3-length sliding window was used, with x i = [x i-2 x i-1 x i ], where the Gaussian kernel was applied. Only the first day, that is 481 samples, was considered for training and estimating the optimal parameters using a 10-fold cross-validation configuration with σ ∈ {2 -5 ; 2 -4 ; • • • ; 2 4 ; 2 5 }. To be comparable with the one-class SVM, we considered |I| = 55 for both methods, obtained from ν = 0.004 in one-class SVM, leading to 60% of common SVs. Table 2 provides a comparative study, with training error given by the cross-validation, and the test error estimated on the next 3 days. This shows the relevance of our approach, where we also illustrate the upper bound on the approximation error derived in Theorem 1. 3

CONCLUSION

In this paper, we investigated a new one-class classification method, by using the coherence criterion. We derived an upper bound on the error of approximating the center of the sphere by the resulting reduced model. We incorporated this criterion into a new kernel-based one-class algorithm by solving a least-squares optimization problem, and considered an incremental and decremental schemes. Experiments were conducted on real datasets to compare our approach to existing methods. Perspectives include the use of this approach to derive an online one-class algorithm for novelty detection.


Jie Chen
email: jie.chen@utt.fr

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

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

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

Paul Honeine
email: paul.honeine@utt.fr

ONLINE SYSTEM IDENTIFICATION UNDER NON-NEGATIVITY AND ℓ 1 -NORM CONSTRAINTS ALGORITHM AND WEIGHT BEHAVIOR ANALYSIS

come

INTRODUCTION

Information processing with ℓ 1 -norm constraint has been a topic of considerable interest during the last five years. Many methods in statistical learning, compressive sensing and related fields have been formulated as convex optimization problems subject to ℓ 1 -norm constraint on the set of parameters to estimate. One motivation for using ℓ 1 -norm constraint is that it produces sparse solutions [START_REF] Candes | An introduction to compressive sensing[END_REF]. Non-negativity constraints are also desired properties that can usually be imposed due to inherent physical characteristics of real-life phenomena. For instance, in the study of a concentration field or a thermal radiation field, non-negativity of observations is a property that can be advantageously exploited by sensor networks [START_REF] Chen | System identification under nonnegativity constraints[END_REF]. Combining non-negativity and ℓ 1 -norm constraints has been addressed within various contexts including, e.g., image deblurring in astrophysics [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF] and hyperspectral data analysis [START_REF] Dobigeon | Bayesian separation of spectral sources under non-negativity and full additivity constraints[END_REF].

In this paper, we consider the problem of minimizing a convex cost function subject to non-negativity constraint and constant ℓ 1 -norm constraint, called hereafter full-additivity constraint, within the context of online system identification. Popular approaches are the projected-gradient type methods and the regularizationbased methods. The former project the results onto the feasible region at each iteration [START_REF] Duchi | Efficient projections onto the l 1 -ball for learning in high dimensions[END_REF][START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF]. The latter minimize the objective function with additional penalty terms that are related to the constraints [START_REF] Rao | Learning mappings in brain machine interfaces with echo state networks[END_REF]. This paper is an extension of our previous work [START_REF] Chen | System identification under nonnegativity constraints[END_REF][START_REF] Chen | Non-negative least-mean-square algorithm[END_REF] in the sense that, in addition to the non-negativity constraint, we also consider now the full-additivity constraint. We propose a strategy that differs from projected-gradient methods and regularization methods in the sense that it integrates the constraints into the coefficient update process. No extra projection procedure or penalty term is required, and constraints are strictly satisfied at each iteration. In order to deal with non-stationary systems, and to reduce computational requirements, we also develop a stochastic gradient algorithm that updates the parameters in an online way. Its mean-weight convergence behavior is analyzed analytically. Simulations show the advantage of the proposed algorithm over some existing methods, and the accuracy of the analysis.

PROBLEM FORMULATION

Consider an unknown system, characterized by a set of real valued discrete-time responses to known stationary inputs. In this paper, we address the problem of deriving a filter ŷ = x ⊤ (n) α(n), subject to constraints (2)-(3), that best models the desired output y(n) corrupted by noise z(n). See Figure 1. Input signal x(n) is zero-mean, and z(n) is a zero-mean additive noise uncorrelated with any other signal. The system identification problem can be formulated as

α o = arg min α J(α) (1)
subject to

α i ≥ 0, i = 1, . . . , L (2)
L i=1

α i = c 0 (3)
with J(α) a convex cost function, L the order of the filter, and c 0 a given positive constant.

WEIGHT UPDATE ALGORITHM

Learning with non-negativity constraints

In [START_REF] Chen | System identification under nonnegativity constraints[END_REF][START_REF] Chen | Non-negative least-mean-square algorithm[END_REF], we have studied an algorithm for system identification involving only the constraint (2), without the full-additivity constraint [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF]. For this problem, using Lagrange multiplier method, the Karush-Kuhn-Tucker conditions at the optimum were combined into the following expression [2, 3]

α o i [-∇ α J(α o )] i = 0 (4)
where the minus sign is just used to make a gradient descent of the criterion J(α) apparent. Equations of

y(n) = x ⊤ (n) α * + z(n) ŷ(n) x(n) Adaptation + - e(n) α(n)
α i (k + 1) = α i (k) + µ i (k) α i (k)[-∇ α J(α(k))] i (5)
where µ i (k) > 0 is a positive step-size that allows us to control convergence. Suppose that α i (k) > 0. The non-negativity of α i (k + 1) is guaranteed if, and only if,

1 + µ i (k) [-∇ α J(α(k))] i > 0. (6)
If [∇ α J(α(k))] i < 0, condition ( 6) is satisfied and the non-negativity constraint does not impose any restriction on the step-size.

Conversely, if [∇ α J(α(k))] i > 0, non-negativity of α i (k + 1) holds if 0 < µ i (k) ≤ 1 [∇ α J(α)] i . (7)
In [START_REF] Chen | Non-negative least-mean-square algorithm[END_REF], we have shown that condition (6) is not necessary for convergence of the algorithm (5) to the solution of problem ( 1)-(2).

Algorithm for fully constrained problem

If non-negativity of α(n) is guaranteed at each iteration, we can make the following variable change to ensure that the full-additivity constraint (3) is satisfied

α j = w j L ℓ=1 w ℓ c 0 . (8)
With w i ≥ 0, the problem becomes unconstrained with respect to constraint (3). The partial derivative of the cost function J with respect to the new variables w i can be expressed as follows

∂J ∂w i = L j=1 ∂J ∂α j × ∂α j ∂w i (9)
with

∂α j ∂w i = c 0 δ ij -α j L ℓ=1 w ℓ (10)
and δ ij the Kronecker symbol. Replacing (10) into ( 9), the partial derivative of J with respect to w i can now be written as

∂J ∂w i = -1 L ℓ=1 w ℓ   -c 0 ∂J ∂α i + L j=1 α j ∂J ∂α j   . (11)
Let us now use the same rule as ( 5) for updating the non-negative entries w i (k). The component-wise update equation is given by

w i (k + 1) = w i (k) + . . . + µ w i (k) L ℓ=1 w ℓ (k)   -c 0 ∂J ∂α i (k) + L j=1 α j (k) ∂J ∂α j (k)   It can be verified that L i=1 w i (k + 1) = L i=1 w i (k)
for all µ by summing both sides from i = 1 to L. This yields

w i (k + 1) = w i (k) + . . . + µ w i (k)   -c 0 ∂J ∂α i (k) + L j=1 α j (k) ∂J ∂α j (k)   .
where 1/ L ℓ=1 w ℓ (k) has been absorbed into µ since it is constant across iterations.

Dividing by L i=1 w i (k + 1) and L i=1 w i (k) the left and right hand sides of the above equation, respectively, and considering the variable change defined by ( 8), we finally obtain

α i (k + 1) = α i (k) + . . . + µα i (k)   -c 0 ∂J ∂α i (k) + L j=1 α j (k) ∂J ∂α j (k)   . (12)
We can verify that

L i=1 α i (k + 1) = L i=1 α i (k)
, which means that the algorithm satisfies the full-additivity constraint as long as the weight vector is initialized by any non-negative vector α(0) whose sum of its components equals c 0 .

Application to mean square error criterion

A general update rule was derived in (12) for any convex cost function J(α). We shall now consider the usual situation where the mean square error is considered, namely

J(α) = E[(α ⊤ x(n) -y(n)) 2 ].
(13) Using criterion (13) in (12) requires that the correlation matrix R x and the intercorrelation vector r xy be estimated. In order to enable online computation, we shall now present a LMS-style algorithm based on the stochastic gradient

-∇J(α) ≈ x(n) e(n) with e(n) = y(n) -x ⊤ (n) α(n).
Substituting this expression into equation (12) yields the following update rule, written in vectorial form,

α(n + 1) = α(n) + . . . + µ diag{α(n)} c 0 x(n)e(n) -1 α ⊤ (n)x(n)e(n) ( 14
)
where 1 is the all-one vector, and diag{•} is a diagonal matrix.

MEAN WEIGHT CONVERGENCE BEHAVIOR

We shall now propose a model to characterize the stochastic behavior of the online algorithm (14). As it is rather challenging to study the performance of such iterative systems, several simplifying assumptions will be made. However, experiments will show that the simulated results and the model outputs match almost perfectly. We define the weight-error vector

v(n) = α(n) -α * (15)
where α * is the solution of the unconstrained problem. See Figure 1. Replacing α(n) by v(n) + α * in (14), and writing e(n) = z(n)v ⊤ (n) x(n), leads us to the following expression

v(n + 1) = v(n) + . . . +µ c 0 z(n) diag{x(n)} v(n) u1 + c 0 z(n) diag{x(n)} α * u2 -c 0 diag{x(n)} v(n) v ⊤ (n) x(n) u3 -c 0 diag{x(n)} α * x ⊤ (n) v(n) u4 -α * α * ⊤ x(n) z(n) u5 -α * v ⊤ (n) x(n) z(n) u6 -v(n) α * ⊤ x(n) z(n) u7 -v(n) v ⊤ (n) x(n) z(n) u8 + v ⊤ (n) x(n) α * α * ⊤ x(n) u9 + v ⊤ (n) x(n) α * v ⊤ (n) x(n) u10 + v ⊤ (n) x(n) v(n) α * ⊤ x(n) u11 + v ⊤ (n) x(n) v(n) v(n) ⊤ x(n) u12
Let us now analyze v(n) in terms of its expected value.

Considering that z(n) is a zero-mean white noise, and using the independence of v(n) and z(n), yields

E[u 1 ] = E[u 2 ] = E[u 5 ] = E[u 6 ] = E[u 7 ] = E[u 8 ] = 0.
After careful calculation, the expected values of the other terms are given by

E[u 3 ] = c 0 diag{R x K(n)} E[u 4 ] = c 0 diag{α * } R x E[v(n)] E[u 9 ] = α * α * ⊤ R x E[v(n)] E[u 10 ] = 1 ⊤ [R x ⊗ K(n)] 1 α * E[u 11 ] = K(n) R x α * E[u 12 ] = K(n) R x E[v(n)]
where

K(n) = E[v(n) v ⊤ (n)]
is the covariance matrix of the weight-error vector, and ⊗ denotes the Hadamard product. It can be noticed that some of the above expected values require second-moments defined by K(n)

in order to update the first-order one E[v(n)]. To simplify the analysis, we use the following separation approximation

K(n) ≈ E[v(n)] E[v ⊤ (n)].
This simplification allows us to study, in an analytical way, the mean-weight behavior of the algorithm. See [START_REF] Hubscher | An improved statistical analysis of the least mean fourth (LMF) adaptive algorithm[END_REF] for a thorough discussion on this approximation. Extensive simulation results have shown us that this approximation provides sufficient accuracy in modeling the mean behavior of the weights. We obtain

E[v(n + 1)] = E[v(n)] + . . . + µ -c 0 diag{R x K(n)} -c 0 diag{α * } R x E[v(n)] + α * α * ⊤ R x E[v(n)] + 1 ⊤ [R x ⊗ K(n)] 1 α * + K(n) R x α * + K(n) R x E[v(n)] (16)
The mean-weight behavior E[α(n)] is then obtained by

E[α(n)] = α * + E[v(n)]. (17)
Simulations will be used to evaluate the performance of this model.

EXPERIMENTAL RESULTS

This section presents simulation examples to illustrate the performance of the proposed algorithm, and to verify the validity of the above analysis. An application to the supervised linear unmixing problem of hyperspectral images is described.

Algorithm behavior illustration

Consider the identification problem in Figure 1, using the least-mean-square criterion. To illustrate the behavior of the proposed algorithm, the impulse response of the system of reference α * , i.e., the solution of the unconstrained problem (1), was set to α * = [0.4 -0.1 0.3 0.05 0 0.01 0.2 -0.01 0.03 0.1] ⊤ Note that α * has two negative entries used to activate the non-negativity constraints, and α * 1 = 1.2. In all the experiments, the same initial vector α(0) > 0 was used, normalized so that the sum of its entries is equal to 1. In the first experiment, the inputs x(n) and the noise z(n) were chosen i.i.d. and drawn from a zero-mean Gaussian distribution with variance σ 2

x = 1 and σ 2 z = 0.01. The non-negativity of the weight coefficients was imposed, and their sum to be equal to 1. The update rule (14) was tested with the step-size values µ = 0.05 and µ = 0.15. The results were averaged over 200 Monte-Carlo runs. Figure 2, first row, shows the learning curves and the mean-weight behavior of the algorithm. One can observe that the algorithm exhibits the very usual behavior of LMS-style approaches, for which small step-size usually means slower convergence but better asymptotic performance. One can also notice that the simulated and theoretical curves representing the evolution of E{α(n)} are perfectly superimposed. Finally, one can check that the non-negativity and the full-additivity constraints are satisfied.

In the second experiment, our algorithm was tested with correlated inputs x(n). We considered an AR model given by x(n) = ax(n -1) + w(n), with a = 0.5. The noise w(n) was chosen i.i.d and zero-mean Gaussian with variance σ 2 w = 1 -0.5 2 , so that σ 2 x = 1 as the first example. The other parameters remained unchanged. As shown in Figure 2, second row, the same conclusions as above can be repeated.

Comparison with another method

The performance of the proposed algorithm was also analyzed in the case where α * is sparse, and compared with the popular projected-gradient scheme. The same experimental setup as above was considered, except that where O M is an all-zero vector of length M . Vector α * was thus containing 100 entries, among which 90 were null. The projected-gradient method described in [START_REF] Duchi | Efficient projections onto the l 1 -ball for learning in high dimensions[END_REF] was considered for comparison purposes. The authors solve the problem (1)-(3) using the update rule

α(n + 1) = Π(α(n) -η∇J(α(n)))
where Π is the Euclidean projection onto the nonnegative ℓ 1 ball. Figure 3(a) compares the learning curves of these two online methods for uncorrelated inputs. Clearly, the proposed method outperformed the projected-gradient method, both in convergence speed and steady-state error. Figures 3(b)-(c) represent the coefficient vector α(n) estimated by the two methods. One can notice that the zeros in α * were more accurately recovered by our algorithm, because zero is a fixed-point for the update equation ( 14), at considerably lower computational complexity. Similar results were obtained with correlated inputs.

Application to hyperspectral data

Our algorithm was used for supervised hyperspectral data unmixing. A key aspect in this problem is to determine abundances (fractions) of different materials under non-negativity and sum-to-one constraints. The studied image is the scene over Moffett Field (CA, USA), captured by the airborne visible infrared imaging spectrometer (AVIRIS). A sub-image of size 50 × 50 pixels was chosen to evaluate the proposed algorithm. This scene is mainly composed of water, vegetation and soil. The spectral signatures of constituent materials were extracted by the VCA algorithm. The linear unmixing problem can be formulated as

α o = arg min α (r -M α) ⊤ (r -M α) subject to α i ≥ 0, i = 1, . . . , L L i=1 α i = 1
where α is the vector of abundances, M the (P × L) matrix of column-wise spectral signatures of constituent materials, and r the observed spectral signature to be unmixed. For this problem, using algorithm (12), we obtained the estimated abundance maps presented in Figure 4. Note that several areas with dominant constituent materials were clearly recovered. These maps are very close to those that would be obtained with active set techniques [START_REF] Lawson | Solving Least Squares Problems[END_REF] or Bayesian methods [START_REF] Dobigeon | Bayesian separation of spectral sources under non-negativity and full additivity constraints[END_REF], which are of a higher level complexity.

CONCLUSION

In this paper, we proposed a method for system identification subject to non-negativity and full-additivity constraints. It reaches an optimal solution by traversing the interior of the feasible region, without any extra projection. It achieved faster convergence and lower steady-state error than a projected-gradient based approach, with significantly less computational complexity. We also studied the mean-weight convergence behavior. Future work will focus on analyzing the convergence of the algorithm in the least-squares sense, and obtaining conditions for convergence.

Figure 1 :

 1 Figure 1: Adaptive system identification under constraints

Figure 2 :

 2 Figure 2: First and second row: results for uncorrelated and correlated input signals, respectively. Left: learning curves; middle and right: Estimated coefficients αi(n) for µ = 0.05 and µ = 0.15, respectively. The red curves represent the model output (17), and the blue curves are averaged Monte-Carlo runs.

α

  * = [0.4 -0.1 O 20 0.3 0.05 0 O 30 0.01 0.2 . . . . . . -0.01 O 40 0.03 0.1] ⊤

Figure 3 :

 3 Figure 3: Comparison between the proposed method and the projected-gradient method. (a) Learning curves. (b) Estimated coefficients αi(n) by the proposed method. (c) Estimated coefficients by the projected-gradient method.

Figure 4 :

 4 Figure 4: Abundance maps estimated by the proposed algorithm (vegetation, soil and water, respectively)


Maya Kallas

Paul Honeine

Cédric Richard

Clovis Francis

Hassan Amoud

Kernel-Based Autoregressive Modeling with a Pre-Image Technique

Keywords: pre-image, kernel machine, autoregressive modeling, pattern recognition, prediction

published or not. The documents may come

INTRODUCTION

Many -if not most -real-life systems are nonlinear by nature. While linear concepts can be easily tackled using simple linear algebra, they fail to adequately explain nonlinear behavior. This is the case of the autoregressive (AR) modeling for time series analysis, where each sample is given by a linear combination of a small number of previous samples. Under the assumption of a (linear) AR process, it is easy to estimate the model parameters, i.e., the weights in the linear expansion, and thus predict future observations from previous ones.

One way to derive nonlinear techniques based on linear ones, is to transform the data with some nonlinear map, and apply the linear algorithm on the transformed data. This is the essence of the kernel-based machines, contributing to the proliferation of nonlinear techniques since Vapnik's Support Vector Machines (SVM) [START_REF] Vapnik | Statistical Learning Theory[END_REF]. The key idea, known as the kernel trick, lies in writing a classical linear algorithm in terms This work was partly supported by the Franco-Lebanese CEDRE program.

of inner products of the transformed data, only to evaluate them using a (positive semi-definite) kernel, without any explicit knowledge of the mapping function. By substituting the inner product with a kernel, the data are implicitly mapped into some high dimension (even infinite-dimension for some kernels) feature space, with essentially no further computational cost. Many nonlinear techniques have been derived on this concept, such as the kernel principal component analysis, kernel Fisher discriminant analysis, and SVM novelty detection, only to name a few.

Mapping the data to the feature space is of great importance to derive nonlinear techniques based on linear ones. Nonetheless, mapping back from the feature space to the input space is also of primary interest. This is mainly because one often needs to interpret the results in the input space, i.e., the signal space in signal processing. Unfortunately, it turns out that the inverse map generally does not exist and only a few elements in the feature space have a valid pre-image in the input space. This is the pre-image problem, as one seeks an approximate solution by identifying data in the input space from its counterpart in the high-dimensional feature space. Many techniques have been proposed in the literature, with a fixed-point iterative method [START_REF] Mika | Kernel PCA and de-noising in feature spaces[END_REF], a method based on the multidimensional-scaling approach [START_REF] Kwok | The pre-image problem in kernel methods[END_REF], or a more direct method based on the relationship between inner-products in both spaces [START_REF] Honeine | Solving the pre-image problem in kernel machines: a direct method[END_REF] (for a recent review, see [START_REF] Honeine | Pre-image problem in kernel-based machine learning[END_REF]).

The linear AR model is one of the most successful, flexible, and easy to use models for the analysis of time series. In this paper, we propose to extend theses advantages to the characterization of nonlinear time series. A natural extension of the linear AR modeling to nonlinear models is derived, in the light of machine learning. To this end, samples are mapped into a nonlinear feature space where, by minimizing the prediction error, the parameters of the nonlinear model are easily estimated using only kernel, without the need to exhibit the nonlinear map. Once the model parameters determined, one can apply a prediction scheme to forecast the future. However, the prediction stage still operates in the feature space. In order to get back to the input space, i.e., the space of samples, we derive an appropriate pre-image technique.

A few attempts have been made to tackle the nonlinear AR model in the light of machine learning literature. A related work by Kumar et al. [START_REF] Kumar | Kernel approach to autoregressive modeling[END_REF] proposes an AR model in the feature space, however, without any ability to predict. Nonlinear modeling and prediction still have not taken full advantage of recent progress in machine learning, although many efforts have been focused to develop nonlinear time series techniques, such as support vector regression [START_REF] Vapnik | Statistical Learning Theory[END_REF], kernel-based Kalman filter [START_REF] Ralaivola | Time series filtering, smoothing and learning using the kernel kalman filter[END_REF], and online prediction with kernels [START_REF] Richard | Online prediction of time series data with kernels[END_REF].

The rest of the paper is organized as follows: In the next section, we carry out a brief description of the AR model. In Section 3, we describe the kernel-based AR model, with both parameter estimation and prediction scheme. Section 4 illustrates the efficiency of the proposed method on several time series.

LINEAR AUTOREGRESSIVE MODEL

The autoregressive (AR) modeling is a well-known prediction method that has been applied successfully in numerous fields. It is defined by a linear prediction formula where each sample in a time series can be predicted from previous samples. Under the assumption of an AR process of order p, a discrete time series x 1 , x 2 , . . . , x n is defined by the model

x i = p j=1 α p-j+1 x i-j , (1)
up to some additive white noise, where the constants α 1 , α 2 , . . . , α p are the model parameters. In other words, each sample is expressed by a linear combination of the p previous samples. To estimate the model parameters from n available samples of the time series, one often minimizes the mean square prediction error, given as

n i=p+1 x i - p j=1 α p-j+1 x i-j 2 .
By setting to zero the derivatives of the above cost function with respect to each α p-j+1 , for j = 1, 2, . . . , p, we get the optimal parameters 1 . Once the model parameters determined, one can predict a future value from the previous samples by using the model (1).

KERNEL-BASED AUTOREGRESSIVE MODELING WITH A PRE-IMAGE TECHNIQUE

The kernel-based AR modeling proposed in this paper combines, on the one hand an AR model in the feature space with parameters determined thanks to the kernel trick, and on the other hand, a pre-image technique to predict the future back into the input space. This is illustrated in Figure 1. 1 Many methods have been proposed to determine the parameters of the AR model, including the use of the Yule-Walker equations and the forwardbackward scheme. Such studies are beyond the scope of this paper. where we identify the AR process (→). To predict back into the input space, a pre-image problem is solved (→).

Autoregressive model in feature space

Let Φ(•) be nonlinear function, mapping any sample x i to an element Φ(x i ) of some feature space. For a time series with samples x 1 x 2 • • • x n , the corresponding elements in the feature space are Φ(x 1 ) Φ(x 2 ) • • • Φ(x n ). We suppose that these elements satisfy an AR process in the feature space, namely

Φ(x i ) = p j=1 α p-j+1 Φ(x i-j ), (2)
where the predicted future is in the feature space. In matrix form, this AR model can be written as

Φ(x i ) = ϕ i α
where ϕ i contains the p previous samples of Φ(x i ),

ϕ i = [Φ(x i-1 ) Φ(x i-2 ) • • • Φ(x i-p )] and α = [α p α p-1 • • • α 1 ] ⊤
the corresponding vector of parameters.

To estimate the parameters α, we minimize the mean square prediction error in the feature space, between the estimated value p j=1 α p-j+1 Φ(x i-j ) and the real one mapped to Φ(x i ). For a sequence of n available samples, we minimize with respect to α the cost function

ξ(α) = n i=p+1 Φ(x i ) - p j=1 α p-j+1 Φ(x i-j ) 2 (3)
where Φ(x i ) 2 = Φ(x i ) ⊤ Φ(x i ) = κ(x i , x i ). This leads to the following expression

ξ(α) = n i=p+1 α ⊤ ϕ ⊤ i ϕ i α -2α ⊤ ϕ ⊤ i Φ(x i ) + Φ(x i ) ⊤ Φ(x i ) .
By taking the derivative of this expression with respect to α,

∇ α ξ(α) = 2 n i=p+1 α ⊤ ϕ ⊤ i ϕ i -ϕ ⊤ i Φ(x i ) ,
and setting it to zero, we get the optimal parameters α of the nonlinear AR model with

α = n i=p+1 ϕ ⊤ i ϕ i -1 n i=p+1 ϕ ⊤ i Φ(x i )
In matrix form, we define the p-by-p matrix K by taking all the inner products between the p previous elements, ϕ ⊤ i ϕ i , with

K = n i=p+1 ϕ ⊤ i ϕ i ,
and the p-by-1 vector k corresponding to

k = n i=p+1 ϕ ⊤ i Φ(x i ).
This leads to a more compact form for estimating the parameters, as follows

α = K -1 k. (4)
The parameter vector α can be estimated using only inner products between pairs of elements in the feature space, defined by the kernel function κ(x i , x j ) = Φ(x i ) ⊤ Φ(x j ). It is clear that K is the p-by-p matrix whose (j, k)-th entry is n i=p+1 κ(x i-j , x i-k ), and k is the p-by-1 column vector whose j-th entry is n i=p+1 κ(x i-j , x i ). We can therefore consider any off-the-shelf (positive semi-definite) kernel to provide a nonlinear AR model. From the learning machines literature, the most used kernel function is the Gaussian radial basis function, of the form κ(x i , x j ) = exp( -1 2σ 2 x ix j 2 ) where σ is the bandwidth of the kernel.

Prediction scheme using a pre-image technique

Once, we have estimated the model parameters α, the prediction stage consists of predicting a new value x i , for i > n, from the p previous samples, x i-1 , x i-2 , . . . , x i-p . By applying the nonlinear AR model (2), we obtain

ψ i = p j=1 α p-j+1 Φ(x i-j ), (5)
where the prediction ψ i lies in the feature space. However, we are more interested in the predicted sample in the original input space. Thus, we need to map back ψ i from the feature space to the input space.

In general, the exact pre-image may not exist, and even if it exists, it may not be unique. This is referred to as the preimage problem, where one identifies the best x * in the input space whose image Φ(x * ) is as close as possible to ψ i . This optimization problem consists of minimizing the distance between elements in the feature space, namely

x * i = arg min x 1 2 ψ i -Φ(x) 2 .
Many methods have been introduced in literature to solve this nonlinear optimization problem. We propose to use the iterative fixed-point method, in the same spirit as in [START_REF] Kallas | Nonlinear feature extraction using kernel principal component analysis with non-negative preimage[END_REF]. By injecting the model (5) into the above expression, we obtain the following optimization problem for any x i :

x * i = arg min x 1 2 p j=1 α p-j+1 Φ(x i-j ) -Φ(x) 2
This optimization problem can be written as

x * i = arg min x J i (x)
where J i (x) is the cost function defined by

J i (x) = - p j=1 α p-j+1 κ(x i-j , x) + 1 2 κ(x, x).
In this expression, the term independent of x, i.e.,

p j=1 α p-k+1 α p-j+1 κ(x i-k , x i-j ), has been removed.

To solve this problem, one may study the gradient of the cost function J i (x) with respect to x. At the optimum, the gradient with respect to x disappears, namely ∇ x J i (x) = 0. The resulting gradient is given as

∇ x J i (x) = - p j=1 α p-j+1 ∂κ(x i-j , x) ∂x + 1 2 ∂κ(x, x) ∂x . (6)
This is the general form for all kernels. This expression can be further simplified for the wide class of radial kernels, such as the Gaussian kernel. In such cases, κ(x, x) is independent of x, thus ∂κ(x, x)/∂x equals to zero, and only the first term in (6) remains. Therefore, the gradient can be expressed as

∇ x J i (x) = - p j=1 α p-j+1 ∂ exp(-1 2σ 2 x i-j -x 2 ) ∂x = - 1 σ 2 p j=1 α p-j+1 κ(x i-j , x) (x i-j -x).
Setting this gradient to zero at the optimum x * i , we get the fixed-point iterative expression

x * i = p j=1 α p-j+1 κ(x i-j , x * i ) x i-j p j=1 α p-j+1 κ(x i-j , x * i )
.

This result can be interpreted as an AR model, in the same spirit as (1), although the parameters are no longer constants, since we have the form x * i = p j=1 β p-j+1 x i-j , with

β p-j+1 = p i=1 α p-i+1 κ(x i-i , x * i ) -1
α p-j+1 κ(x i-j , x * i ).

EXPERIMENT

In this section, we show the relevance of the proposed method on two real-life time series2 : the Mackey-Glass M G 30 time series, modeling the blood cells production evolution with

dx(t) dt = -0.1 x(t) + 0.2 x(t -τ ) 1 + x(t -τ ) 10
with τ = 30, and the Laser time series from the Santa Fe competition (dataset A). For each time series, the first 300 samples were considered to determine the model parameters, as well as estimating the best value for the order. The mean square prediction error was estimated on the next 300 samples.

To give a well-defined benchmark for comparison, we compare the proposed method to several time series prediction techniques: the linear AR model, the multilayer perceptron with a tanh activation function, the support vector regression [START_REF] Vapnik | Statistical Learning Theory[END_REF] and the nonlinear Kalman filter [START_REF] Ralaivola | Time series filtering, smoothing and learning using the kernel kalman filter[END_REF]. For the kernel machines, the Gaussian kernel has been used, with the proper bandwidth value estimated within the estimation stage. For our approach, the bandwidth is set to σ = 0.3 and σ = 0.015 for the Laser and the M G 30 time series, respectively. The value of order was set to p = 6 for both sequences. Table 1 shows the mean square prediction error, where the values of the methods multilayer perceptron, support vector regression and nonlinear Kalman filter are borrowed from [START_REF] Ralaivola | Time series filtering, smoothing and learning using the kernel kalman filter[END_REF].

CONCLUSION

In this paper, we proposed a kernel-based AR modeling for time series analysis and prediction, inspired by recent advances in machine learning. To this end, we bind, on the one hand an AR model in a feature space with model parameters estimated thanks to the kernel trick, and on the other hand, a prediction scheme back into the input space using a pre-image technique. Experiments on a benchmark of real-life time series confirm the relevance of the proposed method.

This strategy in combining kernel machines and AR models opens the way to the development of a range of diverse nonlinear time series techniques. As future work, we are exploring many possibilities in this direction, such as the use of the Yule-Walker equations, the extension to a vector AR model, as well as a nonlinear ARMA model. The choice of the optimal order remains an open question.

Fig. 1 .

 1 Fig.1. Illustration of the kernel-based AR modeling. Samples are mapped from the input space into the feature space (→), where we identify the AR process (→). To predict back into the input space, a pre-image problem is solved (→).

Table 1 .

 1 The mean square prediction error (MSE) for several time series techniques.

	Laser	M G 30

The M G

and the Laser time series are available from http://www.bme.ogi.edu/∼ericwan/data.html.


Maya Kallas

Paul Honeine

Clovis Francis

Hassan Amoud

NON-NEGATIVE PRE-IMAGE IN MACHINE LEARNING FOR PATTERN RECOGNITION

Many real-life applications are nonlinear by nature. Moreover, in order to have a physical interpretation, some constraints should be incorporated in the signal or image processing technique, such as the non-negativity of the solution. This paper deals with the non-negative pre-image problem in kernel machines, for nonlinear pattern recognition. While kernel machines operate in a feature space, associated to the used kernel function, a pre-image technique is often required to map back features into the input space. We derive a gradient-based algorithm to solve the pre-image problem, and to guarantee the non-negativity of the solution. Its convergence speed is significantly improved due to a weighted stepsize approach. The relevance of the proposed method is demonstrated with experiments on real datasets, where only a couple of iterations are necessary.

INTRODUCTION

Constraints are often required to ensure the physical interpretation of many signal and image processing techniques. In pattern recognition on grayscale images, such as deconvolution, deblurring or denoising applications, the result should be a potential valid image. Commonly stored with 8 bits per sampled pixel, each pixel can have 2 8 values within the interval [0, 255], defined by the weakest intensity (0 for black) and the strongest one (255 for white). As the intensity carries the information, numerous applications give rise to grayscale images with a small number of intense pixels, i.e., many pixels have null values. This is the case in many imaging problems, spanning areas such as biomedical and astrophysics [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with nonnegativity constraints[END_REF]. For this reason, one often includes non-negativity constraint on the solution, which leads to zero-valued pixels. In early studies, the non-negativity constraint was introduced for signal deconvolution by Thomas in [START_REF] Thomas | A positive optimal deconvolution procedure[END_REF] and Prost et al. in [START_REF] Prost | Discrete constrained iterative deconvolution algorithms with optimized rate of convergence[END_REF]. Image deconvolution was studied by Thomas et al. in [START_REF] Thomas | Utilisation des multiplicateurs de lagrange pour la restauration d'image avec contraintes[END_REF], while Snyder et al. introduced image deblurring in [START_REF] Snyder | Deblurring subject to nonnegativity constraints[END_REF]. More recently, multiplicative algorithms for signal restoration were studied by Lantéri et al. in [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with nonnegativity constraints[END_REF]. Most of these researches only focus on linear systems [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF][START_REF] Chen | System identification under non-negativity constraints[END_REF], while many real-life applications exhibit nonlinear behavior.

Within the past decade or so, kernel-based machines have been increasingly used in machine learning [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF]. Achieving high accuracy with low computational cost [START_REF] Vapnik | Statistical Learning Theory[END_REF], they have been highly successful in solving many nonlinear problems in classification, regression, prediction and pattern recognition, only to name a few. They rely on the kernel trick, which transforms a linear algorithm into a nonlinear one as long it can be expressed exclusively in terms of inner products between data. By using a positive semi-definite kernel instead of the inner product, we implicitly map the data from the input space into a feature space using a nonlinear map function. The kernel is called reproducing kernel and the corresponding feature space is the so-called reproducing kernel Hilbert space (RKHS).

Even though mapping to a feature space is important, it is usually more interesting to study patterns in the input space rather than their counterparts in the RKHS. However, both spaces are not in bijection and very few elements of the latter have a pre-image in the former. In general, the exact pre-image may not exist and, if it exists, it might not be unique. This defines the pre-image problem, as one seeks a element of the input space whose image, by the same kernel function, is as close as possible to some element in the feature space. Many techniques have been presented in literature in order to solve such nonlinear optimization problem. Mika et al. introduced this optimization problem and proposed a fixed-point iterative method in [START_REF] Mika | Kernel PCA and de-noising in feature spaces[END_REF]. In [START_REF] Kwok | The pre-image problem in kernel methods[END_REF], a multidimensional-scaling (MDS) method was introduced, while lately, in [START_REF] Honeine | Solving the pre-image problem in kernel machines: a direct method[END_REF], a more direct method using relationship between inner-products was presented. See [START_REF] Honeine | Pre-image problem in kernel-based machine learning[END_REF] for a recent review with several applications in signal processing.

In [START_REF] Kallas | Nonlinear feature extraction using kernel principal component analysis with non-negative preimage[END_REF], we showed that the pre-image can be defined with a non-negative additivity of its contributions, by writing the pre-image as a linear model with the available data. However, the non-negativity of the pre-image has not been studied before. This paper deals with nonlinear pattern recognition under non-negativity constraint on the pre-image. To this end, we show that a gradient descent/ascent scheme can be prescribed for solving the pre-image problem, which is somewhat surprising due to the non-linearity and non-convexity of the optimization problem. By controlling its stepsize at each iteration, we guarantee the non-negativity of the solution. Moreover, by weighting the stepsize by the actual value, we derive an algorithm that converges faster, where the sparsity of the solution is privileged. A fortuitous side-effect of the proposed strategy is its re-? The rest of the paper is organized as follows: In the next section, we introduce kernel-based machines with the kernel PCA and describe the denoising scheme. In Section 3, we study the pre-image problem for pattern recognition, illustrated with a gradient approach, while in Section 4, we consider imposing non-negativity constraint on the pre-image. Section 5 gives experimental results illustrating the efficiency of the proposed method on real handwritten digits from the MNIST datasets.

x x * Φ(x) x 1 x 2 x n Φ(x 1 ) Φ(x 2 ) Φ(x n ) ψ 0 ϕ 1 ϕ 2 ϕ m

KERNEL PCA FOR PATTERN RECOGNITION AND DENOISING

Let IR d be the input space to which we associate the Euclidean dot product x i • x j for any x i , x j ∈ IR d . Let κ : IR × IR → IR d be a positive semi-definite kernel, that is i,j α i α j κ(x i , x j ) ≥ 0 for all α i , α j ∈ IR and x i , x j ∈ IR d . The Moore-Aronszajn theorem [START_REF] Vapnik | Statistical Learning Theory[END_REF] states that every positive semi-definite kernel is associated with a unique reproducing kernel Hilbert space H, and vice-versa. This statement is mathematically expressed by a mapping function, Φ : IR d → H, from the input space into the feature space, such that κ(x i ,

x j ) = Φ(x i ), Φ(x j ) H , for any x i , x j ∈ IR d , where • , • H defines the inner product in H.
Next, we study a particular kernel-based machine for nonlinear pattern recognition: kernel PCA, a nonlinear version of the mostly used Principal Component Analysis (PCA). The main idea in PCA is to extract the most relevant directions, and thus the pertinent subspace, from a collection of available data. These directions correspond to the eigenvectors with the highest eigenvalues of the correlation matrix. The same concept can be applied in the feature space, which give rise to the kernel PCA algorithm.

To this end, each data is mapped into a feature space with Φ : IR d → H, where conventional PCA is applied. Let Φ(x 1 ), Φ(x 2 ), . . . , Φ(x n ) ∈ H denote the mapped data. In practice, one does not require the explicit form of the map, since most computations can be done using the concept of the kernel trick, i.e., with Φ(x i ), Φ(x j ) H , for i, j = 1, 2, . . . , n. By analogy to PCA, each relevant axis ϕ is the eigenvector of λ ϕ = C Φ ϕ, where C Φ is the correlation matrix between the mapped data, namely

C Φ = 1 n n j=1 Φ(x j ) Φ(x j ) ⊤ .
In addition, each eigenvector ϕ lies in the span of the Φ-images, thus can be defined by some coefficients α 1 , α 2 , . . . , α n such that ϕ = n i=1 α i Φ(x i ). It is easy to see that these coefficients are obtained by solving the eigen-problem:

λ α = 1 n K α, where α = [α 1 α 2 • • • α n ] ⊤
and K is a n × n matrix defined by κ(x i , x j ) = Φ(x i ), Φ(x j ) H , for i, j = 1, 2, . . . , n. This allows to construct the most relevant subspace in the feature space H, without the need to exhibit any of its elements.

In many pattern recognition applications, one needs to have an access on some elements of the RKHS. Consider the denoising problem using kernel PCA, where the relevant subspace is assumed to be a denoised subspace. Let x be a data corrupted by noise, and let Φ(x) be its image in the feature space. To provide a denoised version of the latter, it is projected onto the relevant subspace. Thus the resulting projection can be written as a linear expansion in terms of the n Φ-images, namely

ψ = n j=1 γ j Φ(x j ), (1)
where γ j = m k=1 n i=1 α k,i α k,j κ(x, x i ), with only m eigenvectors being retained [START_REF] Schölkopf | Kernel pca pattern reconstruction via approximate pre-images[END_REF]. This is illustrated in Figure 1.

PRE-IMAGE PROBLEM FOR PATTERN RECOGNITION: A GRADIENT APPROACH

Many pattern recognition techniques require, not only the feature (1) in the feature space, but also its counterpart in the input space, e.g., the signal space. Getting back, from the feature space to the input space, is the Table 1: Gradient of the objective function (3) for most commonly used kernels, with respect to x.

Kernel

Expression ∇ x J(x)

Polynomial κ p (x i , x j ) = (c + x i • x j ) p n i=1 γ i p κ p-1 (x i , x) x i + p κ p-1 (x, x) x Sigmoid κ S (x i , x j ) = tanh(c (x i • x j ) + σ) n i=1 γ i (1 -κ 2 S (x i , x)) c x i + c (1 -κ 2 S (x, x)) x Exponential κ E (x i , x j ) = exp( 1 σ (x i • x j )) 1 σ n i=1 γ i κ E (x i , x) x i + 1 σ κ E (x, x) x Gaussian κ G (x i , x j ) = exp(-1 2σ 2 x i -x j 2 ) 1 σ 2 n i=1 γ i κ G (x i , x) (x i -x)
pre-image problem. This principle has shown its relevance in many signal processing techniques; see [START_REF] Honeine | Pre-image problem in kernel-based machine learning[END_REF] for a recent review. More precisely, we seek an element of the input space whose image in the feature space is defined by ψ = n j=1 γ j Φ(x j ). The element is the so-called preimage of ψ. Nonetheless, such pre-image may not exist or it may not be unique. Thus, we seek an approximate pre-image x * whose image Φ(x * ) is as close as possible to ψ. This is the pre-image problem, defined by minimizing the distance between the corresponding features in the RKHS, namely

x * = arg min x ψ -Φ(x) 2 H , (2)
where • H denotes the norm in the RKHS.

In the general case, this optimization problem can be written as

x * = arg max

x J(x),

where J(x) is the objective function defined by

J(x) = n i=1 γ i κ(x i , x) - 1 2 κ(x, x). (3)
and the term 1 2 n i=1 n j=1 γ i γ j κ(x i , x j ) is removed since it is independent of x. At the optimum, x * , the gradient of the aforementioned objective function goes to zero. Let ∇ x J(x) denotes this gradient, namely

∇ x J(x) = n i=1 γ i ∂κ(x i , x) ∂x - 1 2 ∂κ(x, x) ∂x . (4)
This is a general form for all kernels. This expression along with the objective function itself vary with the type of the kernel. Let us first consider the well-known Gaussian kernel. This so-called radial kernel, depending on the Euclidean distance x ix j , is defined by

κ G (x i , x j ) = exp(- 1 2σ 2 x i -x j 2 )
for any x i , x j ∈ IR d , where σ is a positive bandwidth parameter. The gradient, with respect to x, of the resulting objective function is given by

∇ x J(x) = n i=1 γ i ∂ exp(-1 2σ 2 x i -x 2 ) ∂x = 1 σ 2 n i=1 γ i κ G (x i , x) (x i -x).
Another used kernel is the polynomial, a direct generalization of the linear inner product. This projective kernel relies on the inner product x i • x j . It is defined by κ p (x i , x j ) = (c + x i • x j ) p , where p ∈ IN + and some positive c. In this case, we have

∇ x J(x) = n i=1 γ i p κ p-1 (x i , x) x i -p κ p-1 (x, x) x,
where κ p-1 (x i , x j ) = (c + x i • x j ) p-1 . Table 1 summarizes the gradient with respect to x of the most commonly used kernels.

We are now in a position to derive gradient-based algorithms to solve the pre-image problem, subject to the non-negativity constraint.

SOLVING THE PRE-IMAGE PROBLEM UNDER NON-NEGATIVITY CONSTRAINTS

We consider the following constrained optimization problem:

x * = arg max

x J(x) subject to x ≥ 0
In this expression, the inequality defines the nonnegativity of each component of the vector x. Off-line optimization techniques are often computationally expensive and not efficient to solve this nonlinear optimization problem. Iterative techniques are essentially based on a fixed-point approach, however only appropriate for convex objective function (see for instance [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with nonnegativity constraints[END_REF]). However, our objective function ( 3) is non-convex, thus the fixed-point technique is not appropriate.

Noisy images

Fixed-point [START_REF] Mika | Kernel PCA and de-noising in feature spaces[END_REF] (100 iterations)

MDS technique [START_REF] Kwok | The pre-image problem in kernel methods[END_REF] Non-negativity with fixed stepsize (6) (100 iterations)

Non-negativity with weighted stepsize (7) (5 iterations) Figure 2: A set of ten "5"-digit images corrupted by a salt-and-pepper noise of density 0.1 (first row), on which we applied the kernel PCA for data denoising. The pre-image results using the fixed-point iterative algorithm [START_REF] Mika | Kernel PCA and de-noising in feature spaces[END_REF] are illustrated (second row), the MDS technique [START_REF] Kwok | The pre-image problem in kernel methods[END_REF] (third row), the non-negative gradient pre-image (6) (fourth row) and the non-negative pre-image with the iterative schema (7) (last row).

Next, we consider an iterative scheme, updating x(t + 1) from x(t). Let [ • ] ℓ denotes the ℓ-th entry operator, and x ℓ (t) the ℓ-th component of x(t), namely [x(t)] ℓ . Component-wise, the gradient ascent scheme is defined by

x ℓ (t + 1) = x ℓ (t) + η ℓ (t) [∇ x J(x(t))] ℓ ,
where η ℓ (t) is a stepsize factor used to control convergence, and [∇ x J(x(t))] ℓ denotes the ℓ-th entry of the gradient of the objective function [START_REF] Prost | Discrete constrained iterative deconvolution algorithms with optimized rate of convergence[END_REF]. A condition on η ℓ (t) must be satisfied to ensure the non-negativity of x ℓ (t + 1) depending on the sign of the gradient. It is easy to see that, when the gradient is positive, there is no restrictions on the stepsize; when it is negative, the stepsize is upper-bounded by

η ℓ (t) ≤ - x ℓ (t) [∇ x J(x(t))] ℓ .
While one can use a stepsize value for each direction in the gradient ascent algorithm, it is often interesting to have a single stepsize value at a given instance t. We shall define the stepsize η(t) such that

η(t) ≤ min ℓ - x ℓ (t) [∇ x J(x(t))] ℓ . (5)
This allows us to write the updating rule in matrix form, with x(t + 1) = x(t) + η(t) ∇ x J(x(t)).

Next, we propose an approach to converge more rapidly toward null values. To this end, the stepsize η ℓ (t) is weighted by the value of x ℓ , which increases the speed to get towards zero. This leads to the expression

x ℓ (t + 1) = x ℓ (t) + η ℓ (t) x ℓ (t) [∇ x J(x(t))] ℓ .
The factor of convergence is now η ℓ (t) x ℓ (t), which results in a new condition on η ℓ (t) for the non-negativity of the solution. For this purpose, we write the above expression as

x ℓ (t + 1) = x ℓ (t) 1 + η ℓ (t)[∇ x J(x(t))] ℓ ,
and obtain the non-negativity condition on 1 + η ℓ (t) [∇ x J(x(t))] ℓ . While no restriction is required when the gradient is positive, in the other case, when [∇ x J(x(t))] ℓ < 0, the stepsize must be upper-bounded by 1/[∇ x J(x(t))] ℓ . In matrix form, the update rule can be written as [START_REF] Chen | System identification under non-negativity constraints[END_REF] where diag[x(t)] is a diagonal matrix whose entries are the elements of vector x(t). To impose the nonnegativity, we consider a single stepsize value η(t) with

x(t + 1) = x(t) + η(t) diag[x(t)] ∇ x J(x(t)),
η(t) ≤ min ℓ - 1 [∇ x J(x(t))] ℓ . (8)
Next, we study the performance of the gradient-based approach, with either the updating rule (7) (resp. ( 6)), and the non-negativity condition (8) (resp. ( 5)).

EXPERIMENTS

In this section, we illustrate the relevance of the proposed method in an image denoising task: real handwritten digits taken from the MNIST datasets1 . We have chosen the digits "5". Each image is defined by grayscale 28 × 28 pixels, with values normalized between 0 and 1. Thus, each image can be written as a 784-dimensional vector. The images were corrupted by adding a Salt-and-Pepper noise, with 0.1 density. A set of 500 images was used to train the kernel PCA with 50 eigenvectors retained. Another set of 10 images, corrupted by the same noise settings, was used to demonstrate the relevance of this denoising technique.

We compared the proposed method to two stateof-the-art techniques: the MDS approach [START_REF] Kwok | The pre-image problem in kernel methods[END_REF], and the fixed-point iterative method [START_REF] Mika | Kernel PCA and de-noising in feature spaces[END_REF]. The latter is defined by the iterative expression

x(t + 1) = n i=1 γ i κ G (x i , x(t)) x i n i=1 γ i κ G (x i , x(t))
. These methods have been successfully applied on many pattern recognition problems, mainly using the Gaussian kernel. For this reason, we consider the same kernel, the bandwidth was fixed for all algorithms (σ = 8). While the maximum number of iterations was set to 100, only 5 iterations were applied for the weighted stepsize algorithm.

We applied the aforementioned techniques along with the proposed algorithms. In Figure 2, the first row shows the noisy images of the digit "5". Results from the proposed method are given in the last two rows. They should be compared to the results obtained by the fixedpoint iterative method and illustrated in the second row, and to the MDS technique illustrated in the third row. As we can see, the fixed point iterative algorithm failed in such applications. The MDS technique did not succeed in recognizing the patterns, as well as the simple gradient with a fixed stepsize. The proposed method, with the weighted stepsize, has shown to be relevant,, the patterns were recognized and all the handwritten digits identified.

CONCLUSION

In this paper, we studied the non-negativity of the preimage in kernel machines. We showed that it is easy to impose the non-negativity using a simple gradientbased scheme. By considering a weighted stepsize in the iterative algorithm, the convergence speed was significantly improved, with remarkable self-regularization property. While the optimization problem is nonlinear and highly non-convex, the proposed technique gave extremely accurate results, with performance illustrated on handwritten digits taken from the MNIST datasets. In future work, we are interested in a general framework for nonlinear non-negative pattern recognition, including this work as well as our previous work [START_REF] Kallas | Nonlinear feature extraction using kernel principal component analysis with non-negative preimage[END_REF]. In the latter, we considered the non-negativity in the coefficients of the linear model that defines the pre-image.

  This work is partly supported under the French-Lebanese collaboration program CEDRE No. 10 SCI F15/L5.

Figure 1 :

 1 Figure 1: Schematic illustration of the pre-image problem for the denoising problem. A given noisy data x is mapped to Φ(x), then projected into the subspace spanned by the most relevant principal axis ϕ 1 , ϕ 2 , . . . , ϕ m . The denoised pattern ψ is mapped back into the input space, to x * .

The datasets are available from the address http://yann.lecun.com/exdb/mnist/.


Jie Chen

Cédric Richard

Henri Lantéri

Céline Theys

Paul Honeine

A GRADIENT BASED METHOD FOR FULLY CONSTRAINED LEAST-SQUARES UNMIXING OF HYPERSPECTRAL IMAGES

Keywords: Hyperspectral imagery, linear unmixing, estimation under constraints

Linear unmixing of hyperspectral images is a popular approach to determine and quantify materials in sensed images. The linear unmixing problem is challenging because the abundances of materials to estimate have to satisfy nonnegativity and full-additivity constraints. In this paper, we investigate an iterative algorithm that integrates these two requirements into the coefficient update process. The constraints are satisfied at each iteration without using any extra operations such as projections. Moreover, the mean transient behavior of the weights is analyzed analytically, which has never been seen for other algorithms in hyperspectral image unmixing. Simulation results illustrate the effectiveness of the proposed algorithm and the accuracy of the model.

INTRODUCTION

Hyperspectral imagery has been applied to a number of areas, including the environment, land use and agricultural monitoring. It can be used for more accurate and detailed information extraction than other types of remotely sensed data, as images provide rich spectral information to identify and distinguish between spectrally similar materials. This technology has received considerable attention.

In hyperspectral imagery, a pixel is a mixture of spectral components associated with a number of pure materials present in the scene. Although nonlinear mixture analysis has begun to find novel applications, linear spectral mixture analysis is still a widely studied approach to determine and quantify materials in sensed images, due to its simpler physical interpretation. In linear unmixing models, measured pixels can be decomposed as linear combinations of constituent spectra, the so-called endmembers, weighted by corresponding fractions named abundances. A key aspect in determining abundances is how to deal with the physical constraints, which are non-negativity and full-additivity. This problem is usually resolved by minimizing a cost function such as the least square criterion subject to these two constraints. In [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF], the authors derive an iterative algorithm referred to as fully constrained least square algorithm, which adopts the sum-to-one constraint in the signature matrix and then performs the active set nonnegative least square algorithm of [START_REF] Lawson | Solving least squares problems[END_REF]. This method does not satisfy the full-additivity constraint exactly, and the active set method for non-negative least square includes substantial computational complexity. In [START_REF] Theys | Linear unmixing of hyperspectral images using a scaled gradient method[END_REF], another class of algorithm is proposed where the full additivity constraint is ensured by projecting estimates at each iteration.

In this paper, we propose a gradient descent method which does not require any extra projection onto the space of constraints. It contains requirements for non-negativity and fulladditivity into the coefficient update process. These two constraints are always fulfilled. Moreover, we derive an analytical analysis of the mean weight transient behavior of the proposed algorithm. To the best of our knowledge, this is unusual in the hyperspectral imaging community. Simulation results illustrate the effectiveness of the proposed algorithm and the accuracy of the model.

LINEAR MIXTURE MODEL

Linear spectral mixing model is largely used in hyperspectral imagery for determining and quantifying materials in sensed images. The physical assumption underlying this model is that each incident photon interacts with one Earth surface component only, and reflected spectra that are collected are not mixed with one another. Let L be the number of spectral bands, and r = [r 1 , r 2 , . . . , r L ] ⊤ an L-by-1 column pixel of an hyperspectral image. Assume that M is the L-by-R target matrix of endmembers, denoted by

[m 1 , m 2 , . . . , m R ],
where each column m i represents an endmember signature. Let α = [α 1 , α 2 , . . . , α R ] ⊤ be a R-by-1 abundance column vector associated with the pixel r. The linear mixing model is given by r

= M α + n ( 1
)
where n is an additive white noise sequence. In order to estimate the abundance vector α, we consider here the leastsquares cost function

J(α) = (r -M α) ⊤ (r -M α) (2)
To be physically meaningful, this model is subject to two constraints on α. The non-negativity constraint requires all the abundances to be nonnegative, that is, α i ≥ 0 for all i. The full additivity constraint, also referred to as sum-to-one constraint, requires

R i=1 α i = 1.
The linear unmixing problem is cast as the following constrained optimization problem

α * = arg min α J(α) (3)
subject to

α i ≥ 0, i = 1, . . . , R (4)
R i=1 α i = 1. (5)

CONSTRAINED ALGORITHM WITH GRADIENT DESCENT METHOD

Solution with non-negativity constraints

In [START_REF] Chen | System identification under non-negativity constraints[END_REF], we studied an algorithm for system identification involving only constraint (4), without the full-additivity constraint [START_REF] Hubscher | An improved statistical analysis of the least mean fourth (LMF) adaptive algorithm[END_REF]. For this problem, using Lagrange multiplier method, the Karush-Kuhn-Tucker conditions at the optimum were combined into the following expression

α o i [-∇ α J(α)] i = 0 (6)
where the minus sign is just used to make a gradient descent of criterion J(α) apparent. Equations of the form g(u) = 0 can be solved with a fixed-point algorithm, under some conditions on function g, by considering the problem u = u+g(u).

Implementing this fixed-point strategy with (6) leads us to the component-wise gradient descent algorithm

α i (k + 1) = α i (k) + µ i (k) α i (k)[-∇ α J(α(k))] i (7)
where µ i (k) > 0 is a positive step size that allows to control convergence. Suppose that α i (k) > 0. Non-negativity of α i (k + 1) is guaranteed if, and only if, 8) is satisfied and nonnegativity constraint does not impose any restriction on the step size. Conversely, if

1 + µ i (k) [-∇ α J(α(k))] i > 0 (8) If [∇ α J(α(k))] i < 0, condition (
[∇ α J(α(k))] i > 0, non-negativity of α i (k + 1) holds if 0 < µ i (k) ≤ 1 [∇J(α)] i (9)

Algorithm for fully constrained problem

If non-negativity of the parameters is guaranteed at each iteration, we can make the following variable change to ensure that the sum-to-one constraint ( 5) is satisfied

α j = w j R ℓ=1 w ℓ (10)
With w i ≥ 0, the problem becomes unconstrained with respect to constraint [START_REF] Hubscher | An improved statistical analysis of the least mean fourth (LMF) adaptive algorithm[END_REF]. The partial derivative of the cost function J with respect to the new variables w i can be expressed as follows

∂J ∂w i = R j=1 ∂J ∂α j × ∂α j ∂w i (11)
where

∂α j ∂w i = ∂wj ∂wi R ℓ=1 w ℓ - ∂ R ℓ=1 w ℓ ∂wi w j R ℓ=1 w ℓ 2 = δ ij -α j R ℓ=1 w ℓ (12)
The Kronecker symbol δ ij results from the derivative ∂w j /∂w i . Replacing (12) into (11), the negative partial derivative of J with respect to w i can now be written as

- ∂J ∂w i = 1 R ℓ=1 w ℓ   - ∂J ∂α i - R j=1 α j - ∂J ∂α j   (13)
Let us now use the same rule as (7) for updating the nonnegative entries w i (k). The component-wise update equation is given by

w i (k + 1) = w i (k) + . . . + µ w i (k) M ℓ=1 w ℓ (k)   - ∂J ∂α i (k) - R j=1 α j (k) - ∂J ∂α j (k)   (14)
It can be easily found that R i=1 w i (k + 1) = R i=1 w i (k) for all µ. The factor ( M ℓ=1 w ℓ (k)) -1 is thus constant and can be absorbed into µ. This yields

w i (k + 1) = w i (k) + . . . + µw i (k)   - ∂J ∂α i (k) - R j=1 α j (k) - ∂J ∂α j (k)   (15)
Dividing by R i=1 w i (k + 1) and R i=1 w i (k) the left and right sides of (15), respectively, and considering the variable change defined by (10), we obtain

α i (k + 1) = α i (k) + . . . + µα i (k)   - ∂J ∂α i (k) - R j=1 α j (k) - ∂J ∂α j (k)   (16)
We can verify that R i=1 α i (k + 1) = R i=1 α i (k), which means that the algorithm satisfies the sum-to-one constraint as long as the weight vector is initialized by any vector α(0) such that R i=1 α i (0) = 1. The gradient of the least-squares criterion (2) is given by

∇ α J(α) = M ⊤ r -M ⊤ M α. (17)
Considering now the component-wise update equation ( 16), we get

α(k + 1) = α(k) + . . . + µ diag{α(k)} ∇ α J(α) -1∇ α J(α) ⊤ α(k) ( 18
)
where 1 is the all-one vector, and diag{•} a diagonal matrix

MEAN WEIGHT TRANSIENT BEHAVIOR

In this section, we are interested in studying the transient behavior of the iteration governed by (18). As it is well known, it is rather challenging to study the performance of such iterative system. Therefore several simplifying assumption will be adopted in the analysis process. However, experiments will show that the simulated and predicted performance match almost perfectly. Firstly, we rewrite the expression (18) as

α(k + 1) = α(k) + . . . +µ diag{α(k)}(M ⊤ r -M ⊤ M α(k)) +µ diag{α(k)}1(M ⊤ r -M ⊤ M α(k)) ⊤ α(k) (19)
Instead of studying α(k) directly, we introduce the weight error vector

v(k) = α(k) -α * ( 20
)
with α * = arg min J(α) the solution of the unconstrained problem. At each iteration, the residual error can be written as follows

e(k) = r -M α(k) = n -M v(k) (21)
Subtracting α * from both sides of (19) and considering the above expression, the following iterative equation with respect to the weight error vector v(k + 1) is obtained after some mathematical calculation

v(k + 1) =v(k) + µp 1 -µp 2 -µp 3 + µp 4 ( 22
)
where vectors p 1 to p 4 are defined by

p 1 = diag{v(k)}M ⊤ n + diag{α * }M ⊤ n p 2 = diag{v(k)}M ⊤ M v(k) + diag{α * }M ⊤ M v(k) p 3 = α * α * ⊤ M ⊤ n + α * v(k) M ⊤ n + v(k) α * ⊤ M ⊤ n + v(k) v(k) ⊤ M ⊤ n p 4 = α * α * ⊤ M ⊤ M v(k) + α * v ⊤ (k) M ⊤ M v(k) + v(k) α * ⊤ M ⊤ M v(k) + v(k)v ⊤ (k) M ⊤ M v(k)
Let us now analyze v(k) in terms of its expected value. Considering that n is zero-mean white noise, and using the independence of v(k) and n, the expected values of p 1 and p 3 are equal to 0. The expected values of p 2 and p 4 are given by

E[p 2 ] = M ⊤ M diag{K(k)} + diag{α * }M ⊤ M E[v(k)] E[p 4 ] = α * α * ⊤ M ⊤ M E[v(k)] + α * trace{K(k) M ⊤ M } + K(k) M M ⊤ α * + K(k) M ⊤ M E[v(k)]
where

K(k) = E[v(k) v ⊤ (k)]
is the covariance matrix of the weight error vector. Expected values E[p 2 ] and E[p 4 ] require second-order moments defined by

K(n) in order to up- date E[v(k)].
To simplify the analysis, we use the following separation approximation

K(k) ≈ E[v(k)] E[v ⊤ (k)]
This simplification allows us a more detailed study of the mean weight behavior using analytical methods. See [START_REF] Hubscher | An improved statistical analysis of the least mean fourth (LMF) adaptive algorithm[END_REF] for a thorough discussion on this approximation. Extensive simulation results have shown that this approximation achieves adequate accuracy in modeling the mean behavior of the weights. Finally, we obtain the following update equation for the mean weight-error vector

E[v(k + 1)] = E[v(k)] + . . . -µ M ⊤ M diag{E[v(k)] E[v ⊤ (k)]} + diag{α * }M ⊤ M E[v(k)]}) + µ(α * α * ⊤ M ⊤ M E[v(k)] + α * trace{E[v(k)] E[v ⊤ (k)] M ⊤ M } + E[v(k)] E[v ⊤ (k)] M M ⊤ α * + E[v(k)] E[v ⊤ (k)] M ⊤ M E[v(k)] (23)
The behavior of the estimated abundance vector α(k) is then described by α(k) = α * + v(k). Simulations will be used to evaluate the model performance.

COMPUTER SIMULATIONS

Synthetic mixture of real spectra

The proposed learning algorithm was used for hyperspectral data unmixing, generated by linear combination of three pure materials with abundance α = [0.3 0.5 0.2]. 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. See 23) and simulated data match almost perfectly. These results clearly show that the derived model can be used for design purposes and, in particular, to derive a condition on the step-size for convergence.

For comparison, we generated a 50-by-50 hyperspectral image with abundance vectors α ij uniformly generated in the simplex defined by the positivity and sum-to-one constraints. The SNR of this image was set to 20 dB. We run the proposed algorithm, the fully-constrained algorithm described in [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF], and the projected-gradient algorithm presented in [START_REF] Theys | Linear unmixing of hyperspectral images using a scaled gradient method[END_REF]. The root mean square error

RMSE = ij α ij -αij 2 /N R
was used to compare these algorithms, and the following results were obtained: RMSE proposed = 0.0107 RMSE [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF] = 0.0103 (with δ = 0.1) RMSE [START_REF] Heinz | Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery[END_REF] = 0.0120 (with δ = 1) RMSE [START_REF] Theys | Linear unmixing of hyperspectral images using a scaled gradient method[END_REF] = 0.0228 These results show the competitive performance of the presented algorithm compared with state-of-the-art methods.

Simulation on real data

The studied image is the scene over Moffett Field (CA, USA), captured by the airborne visible infrared imaging spectrometer (AVIRIS). A sub-image of size 50 × 50 pixels was chosen to evaluate the proposed algorithm. This scene is mainly composed of water, vegetation and soil. The endmembers were extracted by the VCA algorithm with R = 3. Figure 2 shows the estimated abundance maps. Note that several areas with dominant endmembers are clearly recovered.

CONCLUSION

This paper studied a new iterative algorithm for solving the problem of linear unmixing of hyperspectral images. Its transient mean weight behavior was analytically analyzed. Simulations conducted on synthetic data and real data have validated the algorithm and the proposed model. Future works will include detailed study of convergence, including convergence condition.

Fig- ure 1 ( 3 Fig. 1 .

 131 Fig. 1. Spectra of the three endmembers used for synthesizing data (left). Estimated abundances as a function of the iteration number, with two different step-sizes: µ = 0.1 (middle), µ = 0.05 (right). The red dashed curves represent the model output (23). The blue curves are the result of Monte-Carlo simulations.

Fig. 2 .

 2 Fig. 2. Abundance maps estimated by the proposed algorithm (vegetation, soil and water, respectively)