Maya Kallas

Clovis Francis

Paul Honeine

Hassan Amoud

Cédric Richard

Modeling electrocardiogram using

Keywords: kernel machines, ECG signals, autoregressive model, nonlinear models, pre-image problem

come

I. INTRODUCTION

Heart disease remains the main cause of death nowadays. There are many different forms of heart disease, such as the heart failure, the narrowing or blockage of the coronary arteries. One way to monitor the heart is the electrocardiography. The record, which is called the electrocardiogram (ECG), shows the series of waves that relate to the electrical impulses which occur during each beat of the heart, and can be detected at the surface of the body. Many studies have been made in order to detect and classify the heart diseases of the ECG signals [START_REF] Linh | On-line heart beat recognition using hermite polynomials and neuro-fuzzy network[END_REF], [START_REF] Osowski | Support vector machinebased expert system for reliable heartbeat recognition[END_REF], and [START_REF] Kao | Electrocardiogram analysis with adaptive feature selection and support vector machines[END_REF].

An ECG is a time series corresponding to the electrical activity of the heart. The autoregressive (AR) model is one of the most useful technique for time series analysis [START_REF] Chevalier | Autoregressive models with timedependent log area ratios[END_REF]. Its concept is the prediction of a sample using a combination of some previous ones. Being a linear combination, the model is thus defined by the coefficients, or the weights in the linear expansion. In order to use such model, one has to identify its model parameters that are the coefficients and its order which is the number of the previous samples. A technique used to estimate its coefficients is the Yule-Walker equations. It is based on the evaluation of the autocorrelation function. Still, this linear technique is inappropriate for nonlinear systems such as modeling ECG.

In order to extend the AR to the case of nonlinear systems, we propose to apply the concept of the kernel methods. The main idea behind using the kernels is to map the data from an input space to a feature space using a nonlinear mapping function. In practice, it is not necessary to explicitly determine the mapping function, one can transform a linear algorithm into a nonlinear one by applying the kernel trick [START_REF] Aizerman | Theoretical foundations of the potential function method in pattern recognition learning[END_REF]. Initially introduced with Vapnik's Support Vector Machines (SVM) [START_REF] Vapnik | Statistical Learning Theory[END_REF], the use of kernels proliferated in machine learning due to their low computational cost with high performances. Recently, some kernel-based methods were considered for the analysis and prediction of time series data [START_REF] Richard | Online prediction of time series data with kernels[END_REF], including the SVM for regression and kernel Kalman filter [START_REF] Ralaivola | Time series filtering, smoothing and learning using the kernel kalman filter[END_REF].

In this paper, we propose to combine the concept of the kernel machine with the simplicity of the AR model for modeling real ECG signal. The kernel AR is determined using the Yule-Walker equations in the feature space. Thus, parameters are estimated using the (lagged) expected kernels. To this end, we use a Gaussian kernel with the AR model estimated using the Yule-Walker equations. Eventhough, the kernel AR estimates the samples in some feature space. Therefore, being in a high dimensional space, one cannot interpret the sample. Thus, it is required to map back the data from the feature space to the input space where the sample can be interpreted. This necessitate the use of a pre-image technique to predict the future sample into a space where it can be observed. Many techniques have been presented in literature, in order to solve the pre-image problem. See [START_REF] Honeine | The pre-image problem in kernel-based machine learning[END_REF] for a recent review. The paper's authors have solved the pre-image problem for biomedical signal processing in [START_REF] Kallas | Nonlinear feature extraction using kernel principal component analysis with non-negative pre-image[END_REF], [START_REF] Kallas | Nonnegative pre-image in machine learning for pattern recognition[END_REF], [START_REF] Kallas | Sur le problème de la pré-image en reconnaissance des formes avec contraintes de nonnégativité[END_REF].

The rest of the paper is organized as follows: in the next section, we define the Yule-Walker equations used to estimate the coefficients of the AR model. In section III, we present the principle of the kernel methods in machine learning, and we extend the use of the AR model to the nonlinear systems. A pre-image technique is described in section IV in order to estimate the future sample. Section V illustrates the efficiency of the proposed technique on some real ECG signal taken from the MIT-BIH Normal Sinus Rhythm Database.

II. YULE-WALKER EQUATIONS FOR AR MODEL

A real-valued stochastic process x i is said to be an autoregressive process of order p, denoted by AR(p) if there exists α 1 , . . . , α p ∈ IR, with α p = 0, and a white noise i such that

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

x i x i-1 x i-(p-1) x i-p α 1 α p α p-1
Fig. 1. Illustration of the AR concept, where the x i is defined using the parameters α k and the p previous x i-k .

Therefore, the set x p+1 , x p+2 , . . . , x n defines a time series, described by the AR model of order p. The concept of the AR model is illustrated in figure 1. Several techniques in the literature allow to determine the parameters α k . One of them is the resolution of the Yule-Walker equations. They entail the recursive computation of the autocorrelation function ρ of the AR(p)-process. The recursion ρ ( s) = p j=1 α j ρ(s -j), for lag s > 0 refers to the Yule-Walker equations. Having ρ(-s) = ρ(s), and taking the empirical counterparts r(j) of the autocorrelation ρ(j), we get a more formal equation

r = Rα, where r = [r(1) . . . r(p)] , α = [α 1 . . . α p ] , and R =      1 r(1) . . . r(p -1) r(1) 1 . . . r(p -2) . . . . . . . . . r(p -1) r(p -2) . . . 1      .
When the p × p matrix R is invertible, we can rewrite the formal equation as R -1 r = α.

III. AR MODEL WITH KERNEL MACHINES

A. Kernel methods

The aforementioned techniques are used in the case of linear data. In order to extend their use to the nonlinear case, we propose the kernel machines. Kernel functions are typically viewed as providing an implicit mapping of points into a highdimensional space, with the ability to gain much of the power of that space.

Let X be an input space with the canonical (Euclidean) dot product x i • x j for any x i , x j ∈ X . Let κ : X × X → H be a symmetric and continuous function, i.e., a kernel. A kernel is positive semi-definite if and only if a matrix K with entries κ(x i , x j ) for any finite subset of X is positive semi-definite, that is i,j α i α j κ(x i , x j ) ≥ 0 for all α i , α j ∈ IR and all x i , x j ∈ X . Based on the Moore-Aronszajn theorem [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], any positive semi-definite kernel guarantees the existence of a unique 1 feature space (or reproducing kernel Hilbert space) H, where κ defines an inner product. Let Φ(•) denotes the 1 Unique, up to an isometry. mapping function from the input space X , into the feature space H, then κ(x i , x j ) = Φ(x i ), Φ(x j ) H , for any x i , x j ∈ X , where • , • H denotes the corresponding inner product in H, and • H its norm.

B. Kernel AR model

Using the concept of the kernel machines, and applying it to the data of some time series, each x i of the series is mapped into the feature space using the mapping function Φ(•), yielding to Φ(x i ). So in the case of the AR model, when extending it to the nonlinear systems, and applying the kernel machines, it will be defined in the feature space by

Φ(x i ) = p j=1 α j Φ(x i-j ) + Φ i .
Figure 2 illustrates the concept of the kernel AR model. Let µ be the mean of the mapped time series, namely

µ = IE[Φ(x i )],
where IE[•] is the expectation of Φ(x i ) over all x i 's. Then,

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

p j=1 α j )µ = IE[ Φ i ]
. By taking the centered data (Φ(x i-τ ) -µ), with lag τ > 0, we get

Φ(x i ) -µ, Φ(x i-τ ) -µ H = Φ i -IE[ Φ i ], Φ(x i-τ ) -µ H + p j=1 α j Φ(x i-j ) -µ, Φ(x i-τ ) -µ H . (2)
Using the kernel trick, we replace the dot product in the feature space by some kernel function, we define the centered version of a kernel κ(•, •) with

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

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

r κ = R κ α, xi-1 xi-2 x i-(p-1) xi-p x * i Φ(•) Φ(xi-1) Φ(xi-2) Φ(x i-(p-1) ) Φ(xi-p) αp αp-1 α1 α2 ψi = p j=1 αj Φ(xi-j ) ? X H
Fig. 2. Illustration of the Kernel AR model concept, for which, the data are mapped from X to H, where the linear AR concept is applied on the mapped data. However, once ψ i is evaluated, a mapping back from H to X is needed to interpret x i .

where

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

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

α = R -1 κ r κ .
In practice, the expectations are estimated over a set of n available samples. The centered version of the kernel is determined using

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

IV. PRE-IMAGE TECHNIQUE AS A PREDICTION OF FUTURE

SAMPLES

When the learning stage has been completed with some n available samples, the model parameters are estimated. The AR model in feature space can thus estimate the future sample by

ψ i = p j=1 α j Φ(x i-j ),
starting with i = n + 1. This evaluation of p samples in the feature space gives a sample in such space. Nonetheless, in a high space, this sample can not be interpreted. Thus in order to predict a future sample in the original space, one needs to map back ψ i from the feature space into the input space. Nevertheless, a combination can lead to sample outside the image space, as shown in figure 2, thus, the sample may not have a pre-image, or if it exists it may not be unique. This is the pre-image problem, where one seeks to find an approximate solution in the input space, whose image by the nonlinear function is as close as possible to the sample estimated in the feature space. Several techniques have been proposed to solve this ill-posed problem. This includes a multidimensional scaling technique, a conformal map, and a learning scheme. See [START_REF] Honeine | The pre-image problem in kernel-based machine learning[END_REF] for a formal definition of the preimage problem, and a recent survey of the literature.

As aforementioned, the pre-image problem is solved when we can find an approximate solution x * i whose image Φ(x * ) is as close as possible to ψ i . This is defined by the following optimization problem:

x * i = arg min x 1 2 Φ(x) - p j=1 α j Φ(x i-j ) 2 H .
In general, we can write it as

x * i = arg min x J i (x), with J i (x) = - p j=1 α j κ(x i-j , x) + 1 2 κ(x, x), (4)
where the term independent of x has been removed. One technique used to solve such optimization problem is the fixed-point iterative technique. Based on the gradient Fig. 3. The modeling of 2 signals, using the kernel AR method, and the linear AR method. The blue signal presents the original one, the red is the one modeled using the kernel method, and the cyan is the one modeled using the linear AR.

descent scheme, and the type of kernel used, we can implement an iterative scheme to identify the pre-image [START_REF] Mika | Kernel PCA and de-noising in feature spaces[END_REF]. The Gaussian kernel is the most used one. It is defined by

exp(-x i -x j 2 /2σ 2 ),
for any x i , x j ∈ X , where σ defines its bandwidth parameter. The resulting cost function J i (x) is defined by

- p j=1 α j exp(-x i-j -x 2 /2σ 2 ).
Taking the derivative of the above expression with respect to x, and setting it to zero, we get the fixed-point iterative expression

x * i,t+1 = p j=1 α j exp(-x i-j -x * i,t 2 /2σ 2 ) x i-j p j=1 α j exp(-x i-j -x * i,t
2 /2σ 2 ) .

V. EXPERIMENTATIONS

We applied the proposed method to predict future samples of the ECG signals on some real ones taken from the MIT-BIH Normal Sinus Rhythm Database. The ECG signals are considered stationary, foe they are taken in a 1 minute of time. To this end, the gaussian kernel was used. The experiments were divided into two stages: learning and testing. 150 samples from each signal were used in the learning stage, where the order and the coefficients of the AR model and the bandwidth of the kernel were optimized. The same order for each signal is taken for both the linear AR and the kernel AR models. The following 150 samples were used to test the relevance of the proposed kernel AR model, to the linear AR model where its coefficients are estimated using a linear Yule-Walker equations.

A comparison between the kernel AR method, and the linear AR model has been made. The comparison is evaluated based on the calculation of mean square error estimated on the 150 samples of the test stage. The mean square error (MSE) was estimated with

ε = 1 n 2n i=n+1 x * i -x i 2 ,
where x * i is the predicted value at instant i, and x i is the true value of the ECG signal at the same time. Table I shows the MSE for the kernel AR method and the linear AR method for the ECG signals, while the last column shows the mean of the MSE over the 10 signals. As we can see, the kernel AR model gives better results than the linear AR model. The MSE in the first case is approximately 100 times less than in the linear AR case.

Figure 3 shows two examples of modeling signals for normal cardiac rhythm, on the one hand, using the kernel AR technique estimated with the Yule-Walker equations in the feature space, and on the other, using the linear AR technique evaluated with the Yule-Walker equations. The blue line is the original ECG signal, the red is the one modeled using the proposed kernel AR method, and the cyan is the one modeled using the linear AR. As we can see, the proposed method presents a more appropriate modeling than the linear AR method. The linear AR method fails to model the QRS complex of the heart beat, while the kernel AR method models the complex for each person based on his heart rate.

VI. CONCLUSION

Heart disease is one of the main death causes these days. Yet by analyzing and monitoring the ECG signals, one can predict when an arrythmia may occur. In this paper, we have presented a way to predict the evolution of the electrical activity of the heart using kernel machines. To this end, we have combined the concept of the AR model, which is a linear way to predict sample for time series, with the kernel methods in machine learning, where the data are mapped into a high dimensional feature space. We have defined the kernel AR in feature space using the Yule-Walker equations in order to estimate its coefficients. Once the future samples are estimated using this technique, a pre-image technique is required to map the data back to the input space, where the samples can be interpreted. The pre-image technique used was the fixed-point iterative method. The efficiency of the proposed method is illustrated on some real ECG signals taken from the MIT-BIH Normal Sinus Rhythm Database.

As for future work, we consider the estimation of the AR model using the partial autocorrelation function, or the Bayesian Information Criterion. Also, we intend to extend the use of the kernel machines to the autoregressive moving average or autoregressive integrated moving average model.

TABLE I MSE

 I CALCULATED FOR THE LINEAR AR MODEL AND THE KERNEL AR MODEL USING THE YULE-WALKER EQUATIONS TO MODEL ECG SIGNALS, AND THE MEAN OF THE MSE OVER THE 10 SIGNALS.

	Signal number	1st	2nd	3rd	4th	5th	6th	7th	8th	9th	10th	Mean of MSE
												over the 10 signals
	Kernel AR	0.00056	0.00030	0.00021	0.00027	0.00014	0.00022	0.00018	0.00039	0.00021	0.00059	0.000313
	Linear AR	0.0396	0.0091	0.0631	0.0210	0.0203	0.0268	0.0619	0.221	0.0169	0.0077	0.04874

ACKNOWLEDGMENT

This work is partly supported by the Lebanese University and the French-Lebanese research program CEDRE No. 10 SCI F15/L5.


Zineb Noumir

Paul Honeine

Cédric Richard

ONLINE ONE-CLASS MACHINES BASED ON THE COHERENCE CRITERION

Keywords: support vector machines, kernel methods, one-class classification, online learning, coherence parameter

In this paper, we investigate a novel online one-class classification method. We consider a least-squares optimization problem, where the model complexity is controlled by the coherence criterion as a sparsification rule. This criterion is coupled with a simple updating rule for online learning, which yields a low computational demanding algorithm. Experiments conducted on time series illustrate the relevance of our approach to existing methods.

INTRODUCTION

One-class classification for novelty detection has recently generated a great interest in the machine learning community. When time series are considered, an adaptive scheme is required for online detection. In an online learning scenario, training data are available one sample at a time, as opposed to the batch mode where all the samples are presented to the system at the same time. An online learning is also advantageous when dealing with very large datasets. Several applications of interest in signal processing include audio and speech segmentation [START_REF] Davy | Detection of abrupt spectral changes using support vector machines an application to audio signal segmentation[END_REF] and wireless sensor networks [START_REF] Zhang | Adaptive and online one-class support vector machine-based outlier detection techniques for wireless sensor networks[END_REF].

Several methods have been developed to solve the oneclass problem, the most widely studied being the one-class support vector machines (SVM) [START_REF] Tax | One-class classification[END_REF]. It determines a sphere of minimum volume that encloses all (or most of) the available 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 sparseness of the center. The sparsity property states that the center of the sphere explores only a small fraction of the training samples, known as support vectors (SVs). Quadratic programming techniques are often applied to solve this problem. Such approach is inappropriate for online learning.

Many online learning methods have been proposed for SVM in binary classification problems; see for instance [START_REF] Anlauf | The AdaTron: an adaptive perceptron algorithm[END_REF][START_REF] Cauwenberghs | Incremental and decremental support vector machine learning[END_REF] This work was supported by ANR-08-SECU-013-02 VigiRes'Eau. and references therein. So far, there are few attempts to establish online versions of the one-class SVM. In [START_REF] Davy | An Online Support Vector Machine for Abnormal Events Detection[END_REF], it is argued that the binary classification algorithm in [START_REF] Cauwenberghs | Incremental and decremental support vector machine learning[END_REF] cannot be directly implemented for the one-class problem. In [START_REF] Gómez-Verdejo | Adaptive one-class support vector machine[END_REF], a modified formulation of the one-class SVM is presented, following [START_REF] Kivinen | Online learning with kernels[END_REF] where an exponential window is applied to the data. Still, this technique is based on the slow-varying assumption, and several approximations were applied. In another approach, a novelty detection test can also be defined by comparing two one-class SVMs, one trained on the sliding window before the present instance, and one on the sliding window after it. A similar approach is considered in [START_REF] Zhang | Adaptive and online one-class support vector machine-based outlier detection techniques for wireless sensor networks[END_REF] with an application on wireless sensor networks. Nevertheless, this approach is not an online one-class technique.

In order to derive an online version of the one-class SVM machines, the main difficulty remains in the nature of the constrained quadratic optimization problem. Inspired from the least-squares SVM, a one-class technique is proposed in [START_REF] Choi | Least squares one-class support vector machine[END_REF]. However, as pointed out by the authors, it does not have a decision function, and moreover, it loses the sparseness, thus inappropriate for online detection. In this paper, we revisit the one-class problem by considering a least-squares estimation of the center of the sphere. The sparsity of the solution is controlled online by the coherence parameter, borrowed from the literature on sparse approximation problems with dictionaries [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]. The coherence parameter designates the greatest correlation between elements of a dictionary (or two dictionaries) [START_REF] Donoho | Uncertainty principles and ideal atomic decomposition[END_REF][START_REF] Elad | A generalized uncertainty principle and sparse representation in pairs of bases[END_REF]. In signal processing, the quality of representing a signal using a dictionary is studied in [START_REF] Gilbert | Approximation of functions over redundant dictionaries using coherence[END_REF], while in [START_REF] Gribonval | On the exponential convergence of matching pursuits in quasi-incoherent dictionaries[END_REF] the authors use the coherence parameter for signal processing in multichannel transmissions and for source separation problems. More recently, we have adapted the coherence parameter to a dictionary of kernel functions, and applied it with success for online prediction of time series data with kernel functions [START_REF] Honeine | On-line nonlinear sparse approximation of functions[END_REF][START_REF] Richard | Online prediction of time series data with kernels[END_REF]. The coherence criterion provides an elegant model reduction criterion with a less computationally demanding procedure. In this paper, we associate this criterion with a one-class classification algorithm by solving a least-squares optimization problem.

The rest of the paper is organized as follows. Section 2 outlines the one-class SVM. We present our approach in Sec-tion 3. Algorithms for offline and online learning are described in Section 4. Section 5 provides an experimental study on time series.

ONE-CLASS SVM

Let x 1 , . . . , x n be the set of training samples, et Φ(•) be a nonlinear transformation defined by the use of a kernel, κ(x i , x j ) = Φ(x i ), Φ(x j ) , such as the Gaussian kernel κ(x i , x j ) = exp(x ix j2 /2σ 2 ) where σ is the bandwidth parameter. The one-class SVM 1 finds a sphere, of minimum volume, containing most of the 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
In this expression, ν is 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 KKT optimality conditions, we get c = n i=1 α i Φ(x i ), 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)
It is well known that only a small fraction of the training samples contributes to the above model. These samples, called support vectors (SVs), have non-zero α i 's.

In order to provide an online algorithm for one-class SVM, the main difficulty remains in the constrained optimization problem. For instance, upon arrival of a new sample, the coefficients α i 's should be updated subject to constraints (1). However, the upper bound on the α i 's depends on n, thus should be updated and accordingly the values of these coefficients. This problem is studied in [START_REF] Davy | An Online Support Vector Machine for Abnormal Events Detection[END_REF], where 6 conditions are considered. This increases the computational complexity, which depends on the stability of the solution when new samples are added. 1 There exists another formulation of the one-class classification problem: Let w defines the hyperplane separating the origin from the projections of the samples, in the feature space, then the optimization problem is min r,w,ζ

w 2 + 1 νn n i=1 ζ i + r subject to w, Φ(x i ) ≥ r -ζ i for all i
with the decision function given by w, Φ(xn) . The equivalence between both formulations, i.e., between w and c, is studied in [START_REF] Tax | One-class classification[END_REF].

PROPOSED ONE-CLASS APPROACH

Consider the estimation of the center from n available samples, namely

c n = 1 n n i=1 Φ(x i ). (2)
Then, any new sample x that satisfies Φ(x)c n > r can be considered as an outlier, where the distance is

Φ(x)-c n 2 = 1 n 2 n i,j=1 κ(x i , x j )- 2 n n i=1 κ(x i , x)+κ(x, x).
This formulation is inappropriate for large scale data, and unsuitable for online learning as n grows infinitely. The use of a sparse solution for the center provides a robust formulation, appropriate for simple online and offline learning algorithms. Let c I denotes a sparse model of the center c n by using a small subset of the available samples with

c I = i∈I α i Φ(x i ), (3)
where I ⊂ {1, 2, . . . , n}, and let |I| denotes the cardinality of this subset. The distance of any Φ(x) to c I is

Φ(x)-c 2 = i,j∈I α i α j κ(x i , x j )-2 i∈I α i κ(x i , x)+κ(x, x). (4)
The optimization problem consists of properly identifying this subset and estimating the optimal coefficients α i 's. A joint optimization procedure requires advances techniques, as illustrated above with the one-class SVM. In this paper, we propose a separate optimization scheme:

Step 1 Select the most relevant samples in the expansion (3), by using the "coherence" as a sparsification criterion.

Step 2 Estimate the optimal coefficients α i 's, in the sense of the least-squares sense, namely

min α c n -c I 2 .
Next, we propose two different settings, an offline and an online scheme, and derive appropriate learning algorithms.

OFFLINE AND ONLINE ONE-CLASS METHODS

Coherence parameter

The coherence is a fundamental quantity for characterizing dictionaries in sparse approximation problems [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]. It designates the greatest correlation between the elements of a dictionary. For a dictionary of unit-norm 2 elements, {Φ(x i ) | i ∈ I}, the coherence parameter is defined by [START_REF] Richard | Online prediction of time series data with kernels[END_REF]

µ = max i,j∈I i =j | Φ(x i ), Φ(x j ) | = max i,j∈I i =j |κ(x i , x j )|. (5)
This parameter corresponds to the largest absolute value of the off-diagonal entries in the Gram (kernel) matrix, i.e., matrix with entries κ(x i , x j ) for i, j ∈ I. We say that the dictionary is µ-coherent. Accordingly, an orthonormal basis is 0-coherent, while dictionaries with at least two identical elements are 1-coherent.

In an offline learning scheme, we consider the subset with the least coherence, the number of elements being fixed in advance. One can for instance consider the Gram matrix of all entries, followed by a pruning procedure by removing the entries with the largest off-diagonal values.

Optimal parameters

Once the elements of the expression (3) identified, we estimate the optimal coefficients α i 's for i ∈ I. Let α be the column vector of these coefficients. We consider the minimization of the approximation error c nc I , namely

α = arg min αi i∈I 1 n n i=1 Φ(x i ) - i∈I α i Φ(x i ) 2 .
By taking the derivative with respect to each α k , and setting it to zero, we obtain

1 n n i=1 κ(x k , x i ) = i∈I α i κ(x k , x i ), for each k ∈ I.
Written in matrix form, we get Kα = κ, where K is the Gram (kernel) matrix with entries κ(x i , x j ) for i, j ∈ I and κ is the column vector with entries 1 n n i=1 κ(x k , x i ) for each k ∈ I. The final solution is given by

α = K -1 κ. ( 6
)
This problem is well-posed 3 , since the Gram matrix K is nonsingular. In fact, for a µ-coherent dictionary of |I| elements, the eigenvalues of its Gram matrix are greater than or equal to 1 -(|I| -1)µ. This property was previously derived by two of the authors of this paper [15, Proposition 1].

Online Coherence Criterion

In an online learning scheme, we have a new sample at each time step. The sparsification rule determines whether, at time step n + 1, the candidate Φ(x n+1 ) can be well approximated by a combination of the elements of the dictionary. If not, it is added to the dictionary. The coherence-based sparsification criterion consists of inserting Φ(x n+1 ) in the dictionary provided that its coherence remains below a given threshold µ 0 , namely max

i∈I |κ(x i , x n+1 )| ≤ µ 0 . ( 7
)
3 The well-posedness of our approach show its relevance to other methods that are numerically unstable such as [6, See Section 2.4]. Fig. 1. The online method applied on the "frame" distribution. Samples are given in dots •, and elements of the dictionary are shown with circles •. The contour is given by the distance to the estimated center.

The value of the parameter µ 0 ∈ [0, 1[ determines the level of sparsity. In [START_REF] Richard | Online prediction of time series data with kernels[END_REF]Proposition 2], we have demonstrated that, for a compact subspace, the dimension of the dictionary determined with the above sparsification rule remains finite as n goes to infinity.

Let α n be the coefficients estimated at time step n, and K n and κ n the corresponding (Gram kernel) matrix and vector, respectively. Then, the optimal solution from ( 6) is

α n = K -1 n κ n . (8)

Online update scheme

Upon the arrival of a new sample at time n + 1, we consider the coherence criterion [START_REF] Gómez-Verdejo | Adaptive one-class support vector machine[END_REF], to determine whether the model order remains unchanged or is incremented by including the new element in the dictionary.

First case:

max i∈I |κ(x i , x n+1 )| > µ 0
In this case, the new candidate entry Φ(x n+1 ) is not included in the dictionary. Accordingly, the coefficients α i 's are updated, in order to approximate the discarded new entry by the elements of the the expansion (3). Since the dictionary remains unchanged, we have K n+1 = K n . The only change is within the vector κ n in [START_REF] Kivinen | Online learning with kernels[END_REF], which becomes

κ n+1 = 1 n + 1 n κ n + b
where b is the column vector with entries κ(x i , x n+1 ) for all i ∈ I. We get the updating rule of α n+1 from α n as follows: Second case:

α n+1 = K -1 n+1 κ n+1 = n n + 1 α n + 1 n + 1 K -1 n b.
max i∈I |κ(x i , x n+1 )| ≤ µ 0
In this case, the new candidate Φ(x n+1 ) is included into the dictionary. The number of terms in the model ( 3) is incremented, with the dictionary being determined by I ∪ {n + 1}. This leads to the following expressions of the new Gram matrix

K n+1 = K n b b κ(x n+1 , x n+1 ) .
In order to determine K -1 n+1 , we use the the Woodbury matrix identity to get the inverse of K n+1 from K -1 n :

K -1 n+1 = K -1 n 0 0 0 + 1 c -K -1 n b 1 -b K -1 n 1 ,
where c = κ(x n+1 , x n+1 )b K -1 n b and 0 is a column vector of zeros of appropriate size.

The vector κ n+1 is updated from κ n , with

κ n+1 = 1 n + 1 n κ n + b κ n+1
where κ n+1 = n+1 i=1 κ(x n+1 , x i ). The latter expression requires having all the samples in memory. Still one can overcome this difficulty by considering an instant estimation, with κ n+1 = (n + 1) κ(x n+1 , x n+1 ). It is worth noting that the second case, i.e., the incrementation of the model order, seldom occurs.

Finally, by combining these expressions, we get the following update of α n+1 from α n :

α n+1 = 1 n + 1 nα n + K -1 n b 0 - 1 (n + 1) c K -1 n b 1 [nb α n + b K -1 n b -κ n+1 ].

EXPERIMENTATIONS Toy dataset

In order to illustrate our approach, we considered a 2D toy dataset. The training set consisted of 2 000 samples drawn from a "frame" distribution, as illustrated in Figure 1. For such large scale dataset, classical one-class SVM algorithms cannot be applied. In our experiments, we used the Gaussian kernel, with its bandwidth set to σ = 0.5. We applied the proposed online algorithm, as given in Section 4.4. The coherence threshold4 in ( 7) was set to µ 0 = 0.01, which led to a model with 167 elements (8% of the training data). The distance, computed with (4) and shown with the contours in Figure 1, illustrates the boundary of the (hyper)sphere.

Time series domain description

We 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. Chlorination is the process of adding chlorine to water as a method of water purification to make it fit for human consumption as drinking water. The measurements of chlorine were taken from the public water supplies of the Cannes city in France. We considered 2 days of chlorine concentration measures, with measurements sampled at the rate of a sample every 3 minutes. These time series exhibit large fluctuations due to the variations in water consumption and an inefficient control system. See Figure 2.

To capture the structure of the time series, a 10-length sliding window was used, with

x i = [x i-9 • • • x i-1 x i ],
where the Gaussian kernel was applied and σ = 0.2 for all algorithms. We compared the proposed algorithms with two algorithms from the literature: the one-class SVM algorithm [START_REF] Tax | One-class classification[END_REF] and the adaptive one-class SVM [START_REF] Gómez-Verdejo | Adaptive one-class support vector machine[END_REF]. Figure 2 shows the elements retained in the model, obtained either by pruning as given by the coherence parameter in (5) (offline algorithm) or by using the coherence criterion in (7) (online algorithm). For the online algorithm, the coherence threshold was set to µ 0 = 0.5, which led to a model with 49 elements. To get comparable results, the number of elements was set to 49 for the offline algorithm, and the parameters of one-class SVM and adaptive one-class SVM were set accordingly. It is obvious from Figure 2 that the retained elements (in red) are relevant in describing the fluctuations in the time series. The pertinence of the proposed model update rules are shown with the the decision function (distance computed with (4)), as given in Figure 3. This illustrates that our approach, online and offline, is comparable to the one-class SVM and adaptive oneclass SVM results, with lower computational complexity. See Table 1 for the computational cost of these algorithms. Table 1. Estimated computational time (in seconds) of different algorithms for the real time series.

CONCLUSION

In this paper, we investigated a novel online one-class classification method, by associating the coherence parameter with a least-squares optimization problem. Future works include a convergence study and the use of a step-size parameter.

Fig. 2 .

 2 Fig. 2. The time series, with model elements illustrated (in red), for the proposed offline (left figure) and online (middle figure) algorithms, and the adaptive one-class SVM (right figure).

Fig. 3 .

 3 Fig.3. The decision functions estimated by several algorithms. A threshold can be introduced for novelty detection (not illustrated here).

A unit-norm element satisfies Φ(x) = 1 for every x, namely κ(x, x) = 1; otherwise, substitute κ(x i , x j )/ κ(x i , x i ) κ(x j , x j ) for κ(x i , x j ) in the expression of the coherence parameter.

The coherence threshold plays the same role as the ν parameter in SVM. Preliminary experiments are often conducted to select the appropriate values. A comparative study with different values is not provided due to space restrictions.


Zineb Noumir

Paul Honeine

Cédric Richard

KERNELS FOR TIME SERIES OF EXPONENTIAL DECAY/GROWTH PROCESSES

Keywords: Kernel function, one-class classification, normalization, kernel methods, support vector machines

Many processes exhibit exponential behavior. When kernelbased machines are applied on this type of data, conventional kernels such as the Gaussian kernel are not appropriate. In this paper, we derive kernels adapted to time series of exponential decay or growth processes. We provide a theoretical study of these kernels, including the issue of universality. Experimental results are given on a case study: chlorine decay in water distribution systems.

INTRODUCTION

In a wide variety of situations, many processes exhibit an exponential behavior. This behavior follows from natural processes, when a quantity decreases/increases at a rate proportional to its value. Many quantities depend exponentially on time, either by increasing such as in virus spread, economic growth, nuclear chain reaction, or by decreasing such as given in heat transfer, radioactivity, pharmacology, and chemical reactions, only to name a few. Exponential decay1 also arises in many time series, such as the measure of degradation based on vibration (bearings, dampers, etc), as well as any other characteristic whose evolution represents an evolution to failure [START_REF] Gebraeel | Prognostic degradation models for computing and updating residual life distributions in a time-varying environment[END_REF][START_REF] Kaiser | Predictive maintenance management using sensor-based degradation models[END_REF].

While all these problems are of interest, and without loss of generality, this paper deals with chemical reactions, and more precisely the decay of Chlorine concentration in a water supply network. Chlorine is the most used disinfectant to treat water for microbiological protection in a water supply network [START_REF] Gavrila | Modeling water quality in drinking water distribution systems[END_REF]. It is injected at treatment stations, and chlorine concentrations decay over time essentially as water disinfection takes place, while it must be present at the points of water consumption. Chlorine exhibits an exponential decay behavior [START_REF] Islam | Inverse modeling of chlorine concentration in pipe networks under dynamic condition[END_REF][START_REF] Clark | Predicting chlorine residuals in drinking water: second order model[END_REF] while, at any given node (sensor), it largely fluctu-ates due to non-constant consumer demand, temperature, pH, etc. These fluctuations inhibit the use of an adapted model to characterize the domain of operation. See Figure 1.

By learning the system from available data, machine learning provides an elegant approach to overcome modelbased techniques where a perfect understanding of the system is necessary. This is the idea behind kernel-based machines, pioneered by Vapnik's support vector machines (SVM) for classification and regression [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. Domain description using the so-called one-class SVM is proposed in [START_REF] Tax | One-class classification[END_REF], by solving a quadratic programming problem. These machines share an interesting property: sparsity. In fact, only a small number of the training samples contributes to the final model. These samples, the so-called Support Vectors, describe the diversity within the training data.

The performance of a kernel-based learning machine is highly related to the choice of the kernel function. The kernel determines an implicit (nonlinear) map that transforms the data, from the input space to a feature space, and therefore determines its distribution in the latter space. Conventional kernels include the polynomial, the exponential and the Gaussian kernels. Because there is "no free lunch" in kernel selection [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF], the choice of the kernel should reflect prior knowledge about the problem at hand.

In this paper, we design appropriate kernels for time series of exponential decay (or growth) processes. Conventional kernels, such as the Gaussian kernels, are not appropriate for this class of data. This is illustrated in practice by the large dispersion of these data in the resulting feature space. In consequence, one gets "less regularity" using kernel-based machines. This results in many support vectors, even when a single steady (constant-parameter) process is under investigation. In this paper, we propose a class of kernel functions to overcome these drawbacks. Inspired by the normalized kernel [START_REF] Graf | Classification in a normalized feature space using support vector machines[END_REF], we derive several kernels, and study the generalization ability in the sense of universality, as defined in [START_REF] Steinwart | On the influence of the kernel on the consistency of support vector machines[END_REF].

The rest of the paper is organized as follows. Section 2 presents the time series for exponential decay processes. We derive kernels adapted for these data in Section 3. Section 4 outlines the problem of chlorine concentration in water distribution systems, and provide experimental results on real time series. Conclusion and further directions are given in Section 5. But before, we present the classical one-class SVM for domain description.

One-class SVM for domain description

Let x 1 , x 2 , . . . , x n be a set of samples, and let Φ(•) be a nonlinear transformation defined by some positive semi-definite kernel κ(•, •), i.e., κ(x i , x j ) corresponds to the inner product between Φ(x i ) and Φ(x j ).

The one-class SVM finds a sphere, of minimum volume, containing all (or most of) the samples. Its center c and radius r are obtained by solving the following 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
where, ζ 1 , ζ 2 , . . . , ζ n are non-negative slack variables and ν is 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 the optimal model of the center

c = n i=1 α i Φ(x i ), (1)
where the α i 's are the solution of the so-called dual 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.
The one-class SVM for domain description can also be considered for novelty detection tasks. To this end, any new sample x can be regarded as novel to the domain described by the sphere, namely to the distribution of the training samples, when its distance to the center is beyond some threshold, where the distance is given as

Φ(x)-c 2 = i,j α i α j κ(x i , x j )-2 i α i κ(x i , x)+κ(x, x).
The one-class SVM (derived above), as well as SVM machines for regression and classification [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF], enjoy an important property: sparsity. It is well known that only a small fraction of the training samples contributes to the model [START_REF] Gebraeel | Prognostic degradation models for computing and updating residual life distributions in a time-varying environment[END_REF]. These samples, called Support Vectors (SVs), provide a sufficient representation of the samples, thus describing the diversity of the data. Other algorithms based on different SVs selection criteria have been recently applied with success for the one-class problem in [START_REF] Noumir | One-class machines based on the coherence criterion[END_REF][START_REF] Noumir | One simple one-class classification methods[END_REF]. All these algorithms perform essentially identically, while the performance depends crucially on the choice of the kernel function. There is "no free lunch" in kernel selection [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF]. The choice of the kernel should reflect prior knowledge about the problem at hand. The main contribution of this paper is to design appropriate kernels for time series of exponential decay processes.

PROBLEM STATEMENT

A quantity follows an exponential decay process if it decreases at a rate proportional to its value. Mathematically, it satisfies the differential equation

∂x(t) ∂t = -λ x(t),
where x(t) is the positive quantity at time t and λ is a positive number, often called the decay constant. Equivalently, the solution to this differential equation is

x(t) = x(0) e -λ t , (2)
where x(0) is the initial quantity. When the parameter λ is constant, we can define other parameters, such as the (mean) lifetime 1/λ and the half-life (log 2)/λ of the process.

In real-life applications, the decay parameter λ may be time-varying, often by increasing at some instant. A more general model is

x(t) = x(0) e -λ(t) t , (3)
which allows to incorporate a process varying with time, e.g., having a decay parameter piecewise constant. As opposed to include further models to track the variations, e.g., the nonconstant decay parameter λ(t), we consider in this paper a machine learning approach.

We unfold the time-series into a phase space using timedelay embedding process [START_REF] Packard | Geometry from a Time Series[END_REF], by converting the time-series x(t) into a set of vectors

x i = [x(i -m + 1) • • • x(i -1) x(i)] ⊤ , (4)
for i = m, m + 1, . . ., and where m, the embedding dimension, is fixed in advance. 2 With a set of (vector) samples {x m , x m+1 , x m+2 , . . .} at hand, one can consider classical machine learning techniques for pattern recognition or regression. Any kernel function can be applied, the most primitive one being the inner product,

κ(x i , x j ) = x ⊤ i x j , (5)
as well as the polynomial kernel, κ p (x i , x j ) = (x ⊤ i x j + c) p , for p > 0 and c ≥ 0. These kernels often perform poorly, as opposed to the Gaussian kernel

κ G (x i , x j ) = e -1 2σ 2 xi-xj 2 ,
where σ is the bandwidth parameter and x ix j the Euclidean distance. The consistency of this kernel is "guaranteed" by the concept of universal kernels, as derived in [START_REF] Steinwart | On the influence of the kernel on the consistency of support vector machines[END_REF].

The use of the Euclidean metric, with the conventional inner product or distance, within these classical kernels is not adapted for exponential decay time-series. To illustrate this issue, consider two samples, x i , x j (e.g., two consecutive windows), drawn from the same constant-parameter process. It is easy to see that these measures of similarity and dissimilarity, i.e., x ⊤ i x j and x ix j 2 , highly depend on the initial value, namely another sample x(0). We propose to consider kernels that are independent of the multiplicative scaling x 0 . The functions that satisfy the property of invariance under dilations are also called homogeneous of degree zero.

KERNELS FOR TIME SERIES OF EXPONENTIAL DECAY PROCESSES

A valid (reproducing) kernel function is a function that can be expressed as an inner product in some arbitrary space. By Mercer's theorem, this property is satisfied by positive semidefinite functions. Besides, one can engineer a valid kernel from other kernels by applying the simple rules given in Table 1. Moreover, one can provide valid kernels using a power series (polynomial or a Taylor series expansion) of other kernels, as long as the coefficients are non-negative.

Normalized kernels

As a first natural way to remove the dependence on the initial quantity, i.e., x 0 , we consider the normalized kernel,

κ n (x i , x j ) = κ(x i , x j ) κ(x i , x i ) κ(x j , x j ) , (6)
2 The footnotes in this paper are given only for illustration: We say x i ∼ ǫ λ if entries of x i , as defined in (4), satisfy (2) for some x(0).

Rule Expression 1. Linear combination κ(x

i , x j ) = k γ k κ k (x i , x j ) 2. Positive offset κ(x i , x j ) = κ 0 (x i , x j ) + γ 0 3. Product κ(x i , x j ) = k κ k (x i , x j ) 4. Exponential κ(x i , x j ) = e κ0(xi,xj ) 5. Normalization κ 0 (x i , x j )/ κ 0 (x i , x i ) κ 0 (x j , x j )
Table 1. Some simple rules for engineering a valid kernel from available ones, with

γ k ∈ R + .
where κ(•, •) is the classical inner product defined in [START_REF] Clark | Predicting chlorine residuals in drinking water: second order model[END_REF]. This normalized kernel can be viewed as the (linear) inner product kernel applied after normalizing the data, namely by using x i / x i for every x i . 3The normalization removes the scaling factor by projecting the data onto the unit-radius sphere, then applying the conventional inner-product kernel. This projection operation reduces the dimension of the data, while the classical innerproduct kernel performs poorly in general. Next, we take advantage of the universality of the Gaussian kernel in order to provide more interesting kernels.

Kernels with logarithmic distance

We consider another way to derive kernels independent of the initial quantity, by modifying the metric and using the logarithmic distance instead, with

log x i -log x j ,
where the logarithmic function log(•) is applied to each entry of its (vector) argument. It is worth noting that this expression is independent of the base of the logarithm.

Therefore, one can easily replace the classical Euclidean distance with the above logarithmic distance, which allows us to refashion translation-invariant kernels into scale-invariant ones. This is the case of the Radial Basis Functions, with kernels of the form f ( x ix j 2 ) such as the Gaussian kernel. The Gaussian kernel with logarithmic distance is defined by

κ G,log (x i , x j ) = κ G (log(x i ), log(x j )) = e -1 2σ 2 log xi-log xj 2 . ( 7
)
This kernel can be equivalently obtained by considering the inner product (5) on the log-transformed data, namely κ(log(x i ), log(x j )), and apply the Rule 4 (Exponential) and the Rule 5 (Normalization) successively. This leads to a scaleinvariant Gaussian kernel. 4 The use of a logarithmic scale may sometimes be preferred on the linear scale. This is due to the property that this nonlinear transformation stretches low values, and shorten greater ones. Considering a distribution of samples x 1 , x 2 , . . ., then in the logarithmic scale, samples closer to the origin get further away from the "center" of the distribution 5 . Therefore, these samples are more likely to be outliers when applying the one-class SVM algorithm. This yields an interesting single-side detection property which, for the exponential decay data, corresponds to distinguishing large values of the decay parameter λ, namely when decay goes faster.

Finally, the logarithm can be viewed as a power transformation. By using a positive parameter µ, one can consider the family of power transformations with the so-called Box-Cox transformation, where each x(i) is substituted with:

x(i) µ -1 /µ if µ = 0 log x(i) if µ = 0
This transformation has shown its suitability in many statistical problems. This is the case when one needs to transform not normal-like distributed data into data that does follow approximately a normal distribution, namely turning skew unimodal distributions into nearly symmetric normal-like distributions. The use of this transformation allows us to provide a tunable transformation, with the tunable parameter µ, and therefore give a new class of kernels. However, this is beyond the scope of this paper.

More kernels

The normalized kernel ( 6) is a linear kernel applied on the normalized data. Based on this kernel, we propose two different kernels and study some of their properties.

Definition 1. The logarithmic (normalized) kernel is defined as follows:

κ log,n (x i , x j ) = -log(1 -κ n (x i , x j ))
Definition 2. The inverse hyperbolic tangent (artanh) kernel is defined as follows:

κ artanh,n (x i , x j ) = artanh(κ n (x i , x j ))
4 For x i ∼ ǫ λ i and x j ∼ ǫ λ j , i.e., λ i (resp. λ j ) is constant for all the entries of x i (resp. λ j ), we have

κ G,log (x i , x j ) = e -1 2σ 2 τ |λ i (i-τ )-λ j (j-τ )| 2 .
For the same process, i.e., x i , x j ∼ ǫ λ , we get

κ G,log (x i , x j ) = e -1 2σ 2 λm|i-j| 2
. 5 To show this, consider three values (m = 1 for illustration) with values x(i) < x(j) < x(k), and x(i)x(j) = x(j)x(k) . Then, we get in the logarithmic scale: log(x(i))log(x(j)) > log(x(j))log(x(k)) .

Proposition 3. The logarithmic (normalized) kernel and the artanh kernel are valid reproducing kernels.

Proof. To prove this, we consider the expansion of each function into an infinite series:

-log(1 -κ n (x i , x j )) = ∞ k=0 κ n (x i , x j ) k k , (8)
and

artanh(κ n (x i , x j )) = ∞ k=0 κ n (x i , x j ) 2k+1 2k + 1 , (9)
where we have -1 < κ n (x i , x j ) < 1. The validity of these kernels results from applying rules of Table 1, namely Rule 1 and Rule 3. It is worth noting that these rules were also considered to show that the exponential kernel is a valid kernel, by writing it as an infinite series.

These kernels are illustrated in Figure 2, as a function of the normalized kernel κ n (•, •). The artanh kernel is related to the logarithmic (normalized) kernel, since

artanh(z) = 1 2 (log(1 + z) -log(1 -z))
for any z, and more specifically for z = κ n (x i , x j ). Beyond these simple relations, the logarithmic (normalized) kernel defines a feature space that allows to approximate all continuous functions. This property is studied next, with the concept of "universal kernels".

Universality

The concept of universal kernels is proposed in order to study the generalization ability of machine learning classifiers, independent of the learning scheme. In [START_REF] Steinwart | On the influence of the kernel on the consistency of support vector machines[END_REF], the authors establish the existence of the class of "universal" kernels, as the set of kernels that are consistent for a large variety of classification problems, provided a suitably chosen regularization. Examples of universal kernels include the exponential and Gaussian kernels. This is formalized here for the proposed kernel.

Proposition 4. The logarithmic (normalized) kernel is a universal kernel.

Proof. From the expansion of the logarithmic (normalized) kernel into an infinite series, as given in [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF], this kernel takes the form

∞ k=0 a k (κ n (x i , x j )) k ,
with a k > 0, for all k ≥ 0. Due to Corollary 10 in [START_REF] Steinwart | On the influence of the kernel on the consistency of support vector machines[END_REF] (given here 6 for completeness of this paper), we get the universality of this kernel. 6 Corollary 10 in [START_REF] Steinwart | On the influence of the kernel on the consistency of support vector machines[END_REF]: Let 0 < r ≤ ∞ and f : (-r, r) → lR be a C ∞ -function that can be expanded into its Taylor series in 0, namely f (z) = ∞ n=0 an z n for all -r ≤ z ≤ r. If we have an > 0 for all n, then κ(x i , x j ) = f ( x i , x j ) defines a universal kernel on every compact subset of {x ∈ lR d : x 2 2 ≤ r}.

-2.5

2.5 5 ↑ κ n (x i , x j ) ←-κ artanh,n (x i , x j ) κ log,n (x i , x j ) -→ x ⊤ i xj xi xj 1 1 -1 0.5 5
Fig. 2. Illustrations of some kernels, as a function of

x ⊤ i xj xi xj .

EXPERIMENTATIONS

Chlorine in water distribution systems

A drinking water system is a collection of pipes connecting different nodes, reservoirs, pumps and valves used to transport drinking water to consumers. To regulate the quality of drinking water, all pathogenic microorganisms must be destroyed using a chemical disinfectant. The best widely used and inexpensive disinfectant is chlorine. For this reason, the Environmental Protection Agency (USA) and many countries such as France adopt the injection of chlorine into water supply networks, and enforces regulations regarding the minimum amount of chlorine within a drinking water distribution system. Chlorine is injected at treatment stations, and chlorine concentrations decay over time essentially as water disinfection takes place (and reaction occurring at the pipe walls), while it must be present at the points of water consumption. Chlorine sensors are expensive, costly to implement, and unfeasible to equip all the nodes for observability issues. Thus, control algorithms remain inefficient, making automated closed-loop control not yet widely accepted. Currently, open-loop manual control schemes are often investigated. Treatment stations inject just enough chlorine to maintain a minimum predicted chlorine residual everywhere in the network (typically the allowable minimum is 0.2 mg/l). Occasional offline measurements of chlorine concentration taken within some given nodes are sometimes used to adjust the chlorine input.

It is often assumed that the chlorine distribution within a single pipe is described by ∂x(y, t) ∂t + ν ∂x(y, t) ∂y = -λ x(y, t)

where x(y, t) is the chlorine concentration at position y within the pipe and time t, ν is the water velocity, and λ is the kinetics of water, also called first-order reaction constant. One also assumes that chlorine completely mixes at pipe junctions. For a single pipe, the exact solution of the above differential equation is given by x(y, t) = x(y -νt, 0) e -λ t , when y > νt; otherwise x(y, t) = x(0, t -y/ν) e -λ y/ν . In practice, the problem is more complicated. On the one hand, the first-order reaction parameter λ is a function of several parameters such as the wall reaction rate (which depends on the pipe material, roughness, and age, and cannot be directly measured), and the bulk reaction rate (itself a function of temperature, pH, ...). On the other hand, the water velocity ν is never constant, due to the ever fluctuating consumer demand (varies instantly, within the week, within the seasons, ...), while it is most of the time unknown (nominal values of aggregate water demand are often assumed). Combining all these uncertainties and variabilities, a model-based approach is unrealistic.

Experimental results

In this paper, we study chlorine variation at a single location (sensor), and propose to use domain description using oneclass SVM. The measurements of chlorine were taken from the public water supplies of the city of Cannes, in France. We considered 3 days of chlorine concentration measures, with sampling at the rate of a sample every 3 minutes. This time series exhibits large fluctuations due to the variations in water consumption and an inefficient control system (no closedloop control, only one sensor for chlorine observability). See Figure 1.

To capture the structure of the time series, a sliding window of length m = 10 was used, with Fig. 3. The distance profile of all entries to the center of the one-class. The center is determined by classical one-class SVM algorithm, using (from left to right) the normalized linear kernel, the gaussian kernel with logarithmic distance, the logarithmic (normalized) kernel, and the inverse hyperbolic tangent (artanh) kernel.

x i = [x(i -9) • • • x(i -1) x(i)] ⊤ . κ n κ G,log κ log,n κ artanh,n 1
As applications for the one-class SVM method, we considered the four types of kernel as studied in this paper:

• the normalized kernel, in (6);

• the Gaussian kernel with logarithmic distance, in (7);

• the logarithmic (normalized) kernel, in Definition 1;

• and the inverse hyperbolic tangent (artanh) kernel, in Definition 2.

Figure 3 shows the profile of the distance of data to the center of the one-class, as estimated by the one-class SVM algorithm, where different kernels were used. To get comparable results for all these kernels, the number of support vectors was set to 145, which corresponds to 10% of the samples. It is easy to see that the normalized kernel (left figure) provides the largest fluctuations, as opposed to the other proposed kernels. This illustration can be viewed as a better representation of the time series with the latter kernels. One may also consider a threshold (no illustrated) on the value of the distance, in order to provide a decision rule for one-class detection.

CONCLUSION

In this paper, we studied time series from exponential decay/growth processes, such as the chlorine concentration in a water supply network. We derived several kernels appropriate for these data, and studied some theoretical properties. Experiments conducted on real chlorine time series showed the performance of these kernels using a one-class SVM algorithm. As for future work, we are studying other theoretical properties, as well as the use of the proposed approach for classification and regression tasks. We are also investigating the use of these kernels when measurements from several nodes of the network are available.

Fig. 1 .

 1 Fig. 1. The time series of chlorine variation for three days for a given sensor in the water distribution network of the city of Cannes, in France.

Exponential decay includes exponential growth, when one measures the negative quantity.

It is easy to see that, for x i ∼ ǫ λ i and x j ∼ ǫ λ j , we haveκn(x i , x j ) = τ e -λ i (i-τ )-λ j (j-τ ) τ e -2λ i (i-τ ) τ e -2λ j (j-τ )where the summation is on τ = 0, 1, . . . , m -1, which is independent of the initial value x(0).

This work was supported by ANR-08-SECU-013-02 VigiRes'Eau.


Chafic Saidé

Régis Lengellé

Paul Honeine

Cédric Richard

Roger Achkar

DICTIONARY ADAPTATION FOR ONLINE PREDICTION OF TIME SERIES DATA WITH KERNELS

Keywords: Nonlinear adaptive filters, machine learning, nonlinear systems, kernel methods

During the last few years, kernel methods have been very useful to solve nonlinear identification problems. The main drawback of these methods resides in the fact that the number of elements of the kernel development, i.e., the size of the dictionary, increases with the number of input data, making the solution not suitable for online problems especially time series applications. Recently, Richard, Bermudez and Honeine investigated a method where the size of the dictionary is controlled by a coherence criterion. In this paper, we extend this method by adjusting the dictionary elements in order to reduce the residual error and/or the average size of the dictionary. The proposed method is implemented for time series prediction using the kernel-based affine projection algorithm.

INTRODUCTION

Nonlinear models represent a challenge in many practical situations, for which numerous methods have been considered such as neural networks [START_REF] Haykin | Neural Networks and Learning Machines (3rd Edition)[END_REF] and series expansion methods [START_REF] Schetzen | The Volterra and Wiener theories of nonlinear systems[END_REF][START_REF] Wiener | Nonlinear Problems in Random Theory[END_REF]. Function approximation methods based on reproducing kernel Hilbert spaces (RKHS) are of great importance in kernel-based regression methods such as support vector regression [START_REF] Smola | A tutorial on support vector regression[END_REF][START_REF] Suykens | Least squares support vector machines[END_REF]. Computational requirements for kernel-based methods depend on matrices which size increases with the number of observations. This characteristic makes them unsuitable for online applications. Several methods have been proposed to overcome these calculation costs for online applications. The main idea consists in introducing a new sample to the model if it contributes significantly in reducing the approximation error and, if necessary, removing the element which contributes the least. In [START_REF] Dodd | Function estimation in hilbert space using sequential projections[END_REF][START_REF] Phonphitakchai | Stochastic meta descent in online kernel methods[END_REF], the authors proposed a sparsification rule based on the orthogonal projections, while in [START_REF] Engel | Kernel recursive least squares[END_REF] the approximate linear dependence criterion was considered.

The computational cost is further reduced in [START_REF] Honeine | On-line nonlinear sparse approximation of functions[END_REF][START_REF] Richard | Online prediction of time series data with kernels[END_REF], where the authors introduced the dictionary1 coherence criterion, where the inclusion of a new input data into the dictionary is performed if the dictionary still has a small coherence. They demonstrated that the number of elements in the dictionary remains finite with time. The coherence criterion was coupled with Kernel Affine Projection Algorithm (KAPA) and Kernel Normalized Least Mean Squares (KNLMS) as particular cases.

In all above methods, each element injected into the dictionary remains permanently unchanged, even if the non stationarity makes it later useless for estimating the solution. This is the reason why adaptation of the dictionary appears necessary to obtain a better accuracy and/or a smaller size of the dictionary. In this paper, we study the adjustment of the dictionary elements using an adaptive algorithm to make it more efficient in minimizing the resulting approximation error. We make use of the coherence criterion, which gives us the possibility to reduce, in average, the size of the dictionary. To illustrate the efficiency of the proposed method, we present several experiments using a well known benchmark.

BRIEF REVIEW OF THE MODEL REDUCTION METHOD USING THE COHERENCE CRITERION

Consider an online prediction problem and let un ∈ U the input data vector at time step n and dn ∈ R the corresponding desired output of the model. Let k : U × U → R be a kernel and H is the RKHS associated with it. Due to the representer theorem [START_REF] Kimeldorf | Some results on Tchebycheffian spline functions[END_REF][START_REF] Schölkopf | A generalized representer theorem[END_REF], the predicted model can be written as follows:

ψn(•) = m j=1 αn,j κ(•, uw j ) (1)
where the coefficients αn,j ∈ R and ψn(•) is a real-valued function in the Hilbert space associated with the kernel function κ. The κ(•, uw 1 ), . . . , κ(•, uw m ) form a m-elements subset called the dictionary Dn and m ≪ n is the model order by analogy with linear transverse filters. The response of the model to un at time n is:

ψn(un) = m j=1
αn,j κ(un, uw j )

The main problem in kernel online prediction is the unknown order of the model which increases with time and hence, the necessity to control it at each time step. The coherence criterion is of a great importance to overcome this problem and it is widely used to characterize a dictionary in linear sparse approximation techniques [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF].

In [START_REF] Richard | Online prediction of time series data with kernels[END_REF], the authors proposed to define this parameter for kernelbased models by:

µ = max i =j |κ(uw i , uw j )|
where κ is a unit-norm 2 kernel, i.e., κ(u k , u k ) = 1 for all u k . The parameter µ is the largest absolute value of the off-diagonal entries in the Gram matrix and reflects the largest cross-correlation of the elements of the dictionary. The dictionary is said to be µ-coherent. Note that µ = 0 for an orthogonal basis.

In order to derive a model of the form (1) with the coherence parameter, at each time step n, a candidate function κ(•, un) is introduced into the dictionary if the following condition is satisfied:

max κ(.,uw j )∈Dn |κ(un, uw j )| ≤ µ0 (2)
where µ0 ∈ [0, 1[ is a threshold parameter determining the level of sparsity and the coherence of the dictionary. A very important consequence of the use of the coherence criterion is that the dimension m of the dictionary remains finite as n goes to infinity. Obviously, m increases with µ0.

THE KERNEL AFFINE PROJECTION ALGORITHM

Let Dn be a µ0-coherent dictionary at time step n (Dn satisfies ( 2)), and m be its order. In (1), the optimal solution vector αn = (αn,1

• • • αn,m) t is obtained in accordance with the least-squares problem min αn = dn -Hnαn 2 , ( 3
)
where Hn is a p × m matrix whose (i, j)-th element is κ(un-i+1, uw j ). Let p be the number of inputs/outputs used in (3). This means that, at each time step n, only the p most recent inputs {un, ..., un-p+1} and observations dn = (dn • • • dn-p+1) t are considered [START_REF] Ozeki | An adaptive filtering algorithm using an orthogonal projection to an affine subspace and its properties[END_REF][START_REF] Sayed | Fundamentals of Adaptive Filtering[END_REF][START_REF] Liu | Kernel Adaptive Filtering: A Comprehensive Introduction[END_REF]. When a new input data vector un+1 is fed to the model, one of the following two cases occurs:

• First case: maxj=1,...,m |κ(un+1, uw j )| > µ0 κ(•, un+1
) is not introduced into the dictionary Dn and the solution vector αn+1 is updated as follows:

αn+1 = αn+η H t n+1 (ǫI+Hn+1H t n+1 ) -1 (dn+1-Hn+1αn)
where η is a step-size control parameter and ǫI is a regularization factor.

• Second case: maxj=1,...,m |κ(un+1, uw j )| ≤ µ0 κ(•, un+1) is introduced into the dictionary. Thus, Dn+1 = Dn ∪ {κ(•, uw m+1 )} with uw m+1 = un+1. The solution vector αn+1 is updated according to:

αn+1 = αn 0 + η H t n+1 (ǫI + Hn+1H t n+1 ) -1 × dn+1 -Hn+1 αn 0
These recursions define the Kernel Affine Projection Algorithm (KAPA). See [START_REF] Richard | Online prediction of time series data with kernels[END_REF] for more details.

DICTIONARY ADJUSTMENT

In this section, we describe our approach to adjust the elements of the dictionary Dn to obtain D A n using an adaptive method, by substituting each uw k with an appropriate u A w k , for k = 1...m. The objective is to minimize the quadratic approximation error e 2 n :

en = dn -ψn(un) = dn - m j=1
αn,j κ(un, uw j ).

Since the coherence criterion ( 2) is satisfied at each time instant, before any adaptation, for all pairs of elements in the dictionary, this constraint must remain satisfied after the adjustment, leading to

max i =j |κ(u A w i , u A w j )| ≤ µ0. (4)
The i-th element of the dictionary is adapted using a perturbation in the opposite direction of the gradient of the instantaneous quadratic error with respect to uw i , as follows:

u A w i = uw i -νng w i ∀i = 1...m (5)
where

g w i = -2en αn,i ∇u w i κ(un, uw i ) (6)
is the gradient of the instantaneous quadratic error and νn represents the step size used to adjust all the elements of the dictionary. Then, for any pair of dictionary elements, we obtain

u A w i -u A w j = δu -νnδg ∀i, j = 1...m
where δu = uw iuw j and δg = g w ig w j . Under the coherence constraint, νn cannot be chosen arbitrarily. The problem is to determine an appropriate νn at each time step n in order to adapt the dictionary. The iterative approach to find the best step size is explored in the case of Radial Basis Functions, of the form

κ(ui, uj ) = f ( ui -uj 2 ), (7)
where f ∈ C ∞ . A sufficient condition for this function to be a valid positive-definite kernel is the complete monotonicity of the function

f [17], i.e., ( -1)
k f (k) (r) ≥ 0, ∀r ≥ 0 (8)
where f (k) (•) denotes the k-th derivative of f (•). From ( 7), we get

∇u j κ(ui, uj ) = -2(ui -uj ) f (1) ( ui -uj 2 )
and the coherence condition (4) leads to

f ( δu -νn δg 2 ) ≤ µ0. (9)
It is possible to construct a local model of the kernel function, by approximating it with a Taylor series around νn ∼ 0:

f ( δu-νnδg 2 ) =f ( δu 2 )-2νn(δu t δg-νn δg 2 )f (1) ( δu 2 )+O(νn)
Using this approximation, condition (9) becomes

-2 δg 2 ν 2 n -2νnδu t δg)f (1) ( δu 2 +µ0-f ( δu 2 ) ≥ 0. ( 10
)
The discriminant of this quadratic inequality in νn is:

∆ = δu t δg f (1) ( δu 2 ) 2 + 2 δg 2 f (1) ( δu 2 )(µ0 -f ( δu 2 )
If ∆ < 0, there is no constraint on the value of νn; otherwise, if ∆ ≥ 0, the boundary points of (10) will be νi,jand νi,j+ as follows:

νi,j± = -δu t δg f (1) ( δu 2 ) ± √ ∆ -δg 2 f (1) ( δu 2 )
Since the quadratic expression in (10) must be positive, the interval of possible values for νn is ] -∞, νi,j-] ∪ [νi,j+, +∞[ νi,j-νi,j+

+ | - | +
Obviously, for each pair of dictionary elements (uw i , uw j ), νn = 0 always belongs to the feasible domain, since in this case there is no adaptation of the dictionary, and the last was µ0-coherent. Because the constant term in ( 10) is positive, νi,jand νi,j+ have the same sign 3 . The interpretation of the two bounds for νn is straightforward. When adjusting any two elements of the dictionary, each one according to the gradient of the quadratic error with respect to this element, the bounds (νi,j-, νi,j+) must be satisfied to avoid any overlap of the influence regions of the two considered elements, and thus violation of the coherence constraint (4). (see Fig. 1). We now propose a heuristic for selecting νn. Other heuristics could be considered, but our main objective is to illustrate the efficiency of dictionary adaptation. We initially select a reference step size ν0 > 0, as commonly done for adaptive algorithms with a fixed step size. By considering all the (νi,j-, νi,j+) pairs, νn is selected as follows:

νi,j+ νi,ju

A w i u A w i uw i u A w j -g w i
• if max i,j νi,j+ ≤ 0 ⇒ νn = ν0; • if 0 ≤ min i,j νi,j-≤ ν0 ⇒ νn = min i,j • if 0 ≤ ν0 ≤ min i,j νi,j-⇒ νn = ν0; • if 0 ≤ min i,j (νi,j-) + ≤ ν0 ⇒ νn = min i,j (νi,j-) + ; • if 0 ≤ ν0 ≤ min i,j (νi,j-) + ⇒ νn = ν0.
In these expressions, (νi,j-) + indicates all the positive values of νi,j-'s. Note that ν0 must be selected relatively small. If ν0 is too large, the elements of the dictionary can be spread over a non useful region of the input space, inducing an increase of the size of the dictionary without reducing the approximation error.

EXPERIMENTATIONS

In our experimentation, we used the same benchmark and the same parameter settings as in [START_REF] Richard | Online prediction of time series data with kernels[END_REF]. It consists of a one step prediction of the following discrete-time nonlinear dynamical system:

vn = 1.1 exp(-|vn-1|) + un dn = v 2 n
where un and dn are the input and the desired output, respectively. The data were generated from the initial condition v0 =0.5. The input was sampled from a zero-mean Gaussian distribution with standard deviation 0.25. The system output was corrupted by an additive zero-mean white Gaussian noise with standard deviation equal to 1.

The system to identify is of the form dn = ψn(un).

The KAPA algorithm is used with the following parameters settings: number of Inputs/Outputs p = 3, step-size control parameter η = 0.01, regularization factor ǫ = 0.07 (see section 3 for details).The acronym AKAPA is adopted to indicate the dictionary Adaptation for Kernel Affine Projection Algorithm.

A set of 200 time series of 3000 samples each was used to compare the KAPA and AKAPA using the Normalized Mean Squared Error (NMSE) which was computed over the last 500 samples according to:

NMSE = E 3000 i=2501 (dn -ψn(un)) 2 3000 i=2501 d 2 n
Another indicator that must be computed is the average final size of the dictionary m calculated for the 200 time series. For the simulations, we adopted both the Gaussian kernel κ(ui, uj) = exp(uiuj 2 /2σ 2 ) with a bandwidth σ = 0.42 and Exponential kernel κ(ui, uj ) = exp(uiuj /σ) with a bandwidth σ = 0.33.

The learning curves shown in Figure 2 and Figure 3 depict the Mean Squared Error and compare different settings for the coherence criterion µ0 and for the adjustment step size ν0. ν0 has been selected using a rough grid search so as to obtain the best performances for the given value of µ0. Table 1 gives a summary of the obtained results for the KAPA and the AKAPA algorithms. The obtained results lead to the following observations:

1. For the same coherence parameter µ0 = 0.3, AKAPA shows a 23.92% increase in the average size of the dictionary m with a decrease of 41.60% in the NMSE for the Gaussian kernel (12.12% and 16.93% respectively for the Exponential kernel). See rows 1 and 3 in Table 1. Fig. 3. Learning curves for KAPA and AKAPA using an Exponential kernel (σ = 0.33) for different coherence parameters µ0 and gradient step sizes ν0.

2. Comparing NMSE for the same average sizes of the dictionaries m using different values for µ0, AKAPA led to a decrease of 38.08% for the Gaussian kernel and 15.15% for the Exponential kernel. See rows 2 and 3 in Table 1.

From these observations we can deduce that, if µ0 and ν0 were properly selected, the size of the dictionary and the approximation error can be greatly reduced. These results were also observed on other (synthetic and real) time series, omitted here due to space limitations.

CONCLUSION

In this paper, we demonstrated the interest of adjusting the dictionary elements within the context of online predictions with kernel-based methods. Our idea was to use an iterative gradient adaptation algorithm that satisfies a coherence measure for the elements of the dictionary. Pruning the dictionary to reduce its size, as well as other adaptation algorithms, will be considered in our future work.

Fig. 1 .

 1 Fig. 1. A 2D illustration showing the constraint of choosing νn ≤ νi,jor νn ≥ νi,j+ to avoid the overlap of the influence regions of u A w i and u A w j .

Fig. 2 .

 2 Fig. 2. Learning curves for KAPA and AKAPA using a Gaussian kernel (σ = 0.42) for different coherence parameters µ0 and gradient step sizes ν0.

  AKAPA, µ0= 0.25, ν0 = 0.005 m = 5.41, NMSE = 0.12751 AKAPA, µ0= 0.3, ν0 = 0.01 m = 6.29, NMSE = 0.12786 KAPA, µ0 = 0.3 m= 5.61, NMSE =0.15392 KAPA, µ0 = 0.345 m = 6.26, NMSE = 0.15070

Table 1 .

 1 Experimental Setup and Performance, with p = 3, η = 0.01, and ǫ = 0.07

					Gaussian kernel (σ = 0.42)	Exponential kernel (σ = 0.33)
			Algorithm		Parameter Settings	m	NMSE	Parameter Settings	m	NMSE
			KAPA [10]		µ0=0.3		2.655	0.15970	µ0=0.3	5.61 0.15392
			KAPA [10]		µ0=0.49		3.285	0.14588	µ0=0.345	6.26 0.15070
			AKAPA (this paper)	µ0=0.3, ν0=0.04	3.29	0.090327	µ0=0.3, ν0=0.01	6.29 0.12786
			AKAPA (this paper)	µ0=0.2, ν0=0.033	2.865	0.09401	µ0=0.25, ν0=0.005 5.41 0.12751
	10
	0	500	1000	1500	2000	2500	3000

The term dictionary stands for a set of input vectors (or their corresponding kernel functions in the RKHS) used to estimate the nonlinear model.

Otherwise, replace κ(•, u k ) with κ(•, u k )/ κ(u k , u k ).

This property is due to the fact that the kernel function satisfies the validity condition[START_REF] Engel | Kernel recursive least squares[END_REF], namely f (1) ( δu 2 ) ≤ 0, and µ 0 ≥ f ( δu 2 ) as the coherence constraint is satisfied at any time step.


Zineb Noumir

Paul Honeine

Cédric Richard

ON SIMPLE ONE-CLASS CLASSIFICATION METHODS

The one-class classification has been successfully applied in many communication, signal processing, and machine learning tasks. This problem, as defined by the oneclass SVM approach, consists in identifying a sphere enclosing all (or the most) of the data. The classical strategy to solve the problem considers a simultaneous estimation of both the center and the radius of the sphere. In this paper, we study the impact of separating the estimation problem. It turns out that simple one-class classification methods can be easily derived, by considering a least-squares formulation. The proposed framework allows us to derive some theoretical results, such as an upper bound on the probability of false detection. The relevance of this work is illustrated on well-known datasets.

INTRODUCTION

The one-class classification machines has become a very active research domain in machine learning [START_REF] De Ridder | An experimental comparison of one-class classification methods[END_REF][START_REF] Schölkopf | Estimating the support of a highdimensional distribution[END_REF], providing a detection rule based on recent advances in learning theory. In one-class classification, the problem consists in covering a single target class of samples, represented by a training set, and separate it from any novel sample not belonging to the same class, i.e., an outlier sample. It has been successfully applied in many novelty detection and classification tasks, including communication network performance [START_REF] Zhang | One class support vector machine for anomaly detection in the communication network performance data[END_REF], wireless sensor networks [START_REF] Zhang | Adaptive and online one-class support vector machine-based outlier detection techniques for wireless sensor networks[END_REF], forensic science [START_REF] Ratle | A comparison of one-class classifiers for novelty detection in forensic case data[END_REF], detection of handwritten digits [START_REF] Tax | Kernel whitening for one-class classification[END_REF] and objet recognition [START_REF] Kemmler | One-class classification with gaussian processes[END_REF], only to name a few. Moreover, it has been extended naturally to binary and multiclass classification tasks, by applying a single one-class classifier to each class and subsequently combining the decision rules [START_REF] Hempstalk | Discriminating against new classes: One-class versus multi-class classification[END_REF].

Since only a single class is identified, it is essentially a data domain description or a class density estimation problem, while it provides a novelty detection rule. Different methods to solve the one-class problem have been developed, initiated from the so-called one-class support vector machines (SVM) [START_REF] Tax | One-class classification[END_REF][START_REF] Schölkopf | Estimating the support of a highdimensional distribution[END_REF]. The one-class classification task consists in identifying a sphere of minimum volume that englobes all (or most of) the training data, by estimating jointly its center and This work was supported by ANR-08-SECU-013-02 VigiRes'Eau. its radius. These methods exploit many features from conventional SVM [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF], including a nonlinear extension thanks to the concept of reproducing kernels. They also inherit the robustness to outliers in the training set, by providing a sparse solution of the center. This sparse solution explores a small fraction of the training samples, called support vectors (SVs), and lying outside or on the sphere.

In one-class SVM as defined in [START_REF] Tax | One-class classification[END_REF][START_REF] Schölkopf | Estimating the support of a highdimensional distribution[END_REF], the resulting convex optimization problem is often solved using a quadratic programming technique. Several efforts have been made in order to derive one-class classification machines with low computational complexity [START_REF] Liu | Fast support vector data descriptions for novelty detection[END_REF]. In the same sense as least-squares SVM is derived from the classical SVM method [START_REF] Suykens | Least squares support vector machine classifiers[END_REF][START_REF] Rifkin | Regularized least squares classification[END_REF], some attempts have been made to derive from the one-class SVM a least-squares variant, such as in [START_REF] Choi | Least squares one-class support vector machine[END_REF]. However, unlike the former, the latter do not have a decision function, thus inappropriate for novelty detection.

In this paper, we propose to solve the one-class problem by decoupling the estimation of the center and the radius of the sphere englobing all (or most of) the training samples. In the same spirit as the classical one-class SVM machines, we consider a sparse solution with SVs lying outside or on the sphere. It turns our that the optimal sparse solution can be defined using a least-squares optimization problem, thus leading to a low computational complexity problem. This framework allows us to derive some theoretical results. We give an upper bound on the probability of false detection, i.e., probability that a new sample is outside the sphere defined by the sparse solution.

As opposed to the jointly optimization of the center and radius by the classical one-class SVM approach, our strategy decouples the estimation problem, thus provides a sub-optimal solution. Consequently, the proposed approach should degrade the performance. In practice, we found that the performance is essentially equivalent to the classical technique, while operating a dramatic speed-up. This is illustrated on experiments from a well-known benchmark for one-class machines [START_REF] Liu | Fast support vector data descriptions for novelty detection[END_REF].

The rest of the paper is organized as follows. Section 2 outlines the classical one-class SVM. We describe our approach in Section 3, and derive theoretical results in Section 4. Section 5 illustrates the relevance of our approach on real datasets. Conclusion and further directions are given in Section 6.

CLASSICAL ONE-CLASS SVM

Thanks to the concept of reproducing kernels [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], a (positive semi-definite) kernel function κ(•, •) defines a nonlinear transformation Φ(•) of the input space into some feature space. A sphere defined in the latter corresponds (is preimaged [START_REF] Honeine | Preimage problem in kernel-based machine learning[END_REF]) to a nonlinear characteristics in the input space. It turns our that only the inner product is often required, which can be evaluated using a kernel function, Φ(x i ), Φ(x j ) = κ(x i , x j ) for any x i , x j from the input space X .

The one-class SVM was initially derived in [START_REF] Schölkopf | Estimating the support of a highdimensional distribution[END_REF] for the estimation of the support of a distribution with the ν-SVM, and in [START_REF] Tax | One-class classification[END_REF] for novelty detection with the so-called "support vector data description". The principle idea is to find a sphere, of minimum volume, containing all 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 subject to Φ(x i ) -c 2 ≤ r 2 for i = 1, 2, . . . , n
While the above constraint may be too restrictive, one may tolerate a small fraction of the samples to be outside the sphere. This yields robustness, in the sense that it is less sensitive to the presence of outliers in the training dataset. For this purpose, let ν be a positive parameter that specifies the tradeoff between the sphere volume and the number of outliers. Then the problem becomes the estimation of c, r, and a set of non-negative slack variables

ζ 1 , ζ 2 , . . . , ζ n : 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
By introducing the Karush-Kuhn-Tucker (KKT) 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. (2)
In accordance with the KKT conditions, each sample x i can be classified into three categories: α i = 0 corresponds to a sample lying inside the sphere, samples with 0 < α i < 1 νn lie on the sphere boundary, and samples with α i = 1 νn lie outside the sphere, i.e., are outliers. The samples with non-zero α i are called support vectors (SVs) since they are sufficient to describe the center as defined in expression (4). In practice, only a very small fraction of the data are SV. Let I sv be the set of indices associated to SV, namely

α i = 0 if i ∈ I sv ; α i = 0 otherwise.
Finally, the optimal radius is obtained from any SV lying on the boundary, namely any

x i with 0 < α i < 1 νn , since in this case Φ(x i ) -c = r. This is equivalent to r = min i∈Isv Φ(x i ) -c .
Therefore, the decision rule that any new sample x is not an outlier is given as Φ(x)c < r, where the distance is computed by using

Φ(x)-c 2 = i,j∈Isv α i α j κ(x i , x j )-2 i∈Isv α i κ(x i , x)+κ(x, x).
(3)

SIMPLE ONE-CLASS METHODS

Back to basics, the center (or empirical first moment) of a set of samples is defined by

c n = 1 n n i=1 Φ(x i ), (4)
and the radius of the sphere englobing all (or most of ) the samples can be easily considered, where the distance is evaluated using (3) where α i = 1/n for all i = 1, 2, . . . , n and I sv = {1, 2, . . . , n}. While the sphere defined by the above full-model is extremely sensitive to outliers, one may consider a sparse solution by incorporating a small number of relevant samples in the model. This is essentially the spirit of the classical one-class SVM, which estimates jointly the center and the radius, by identifying the SVs. Our approach towards a sparse solution is based on three steps:

• Determine the full-model center from (4);

• Identify the SVs as the farthest samples from the center;

• Estimate accordingly the sparse model parameters.

Sparsification rule

The classical one-class SVM method provides a sparse model for the center, where only samples outside or lying on the sphere are SVs. Inspired by this result, we consider in our approach the distance criterion to identify this subset.

The set of SVs can be identified by considering the distance of each sample to the center, namely

I = arg max k∈I Φ(x k ) -c n 2 = arg max k∈I -2 n i=1 κ(x i , x k ) + nκ(x k , x k ),
where the number of SVs is fixed in advance. Once the set {x i | i ∈ I} is determined, the radius is given as

r = min i∈I Φ(x i ) -c I ,
where c I is the sparse model of the center defined by

c I = i∈I α i Φ(x i ), (5)
the coefficients α 1 , α 2 , . . . , α n being estimated as follows.

Sparse formulation of the center

Consider the error of approximating c n with the sparse model c I , c nc I , which indicates the wellness of such approximation using a small subset of the training data. The coefficients in ( 5) are estimated by minimizing this error, with α = arg min α1,...,αn

1 n n i=1 Φ(x i ) - i∈I α i Φ(x i ) 2 , (6)
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 , namely -2 Φ(x k ), c n -i∈I α i Φ(x i ) , 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 κ, (7)
where K is the 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. To make this problem well-posed in practice, we include a regularization parameter ν, namely α = (K + νI) -1 κ, where I is the identity matrix of appropriate size. The error of approximating the center with the above solution is

c n -c I 2 = c n 2 -2α ⊤ κ + α ⊤ Kα = 1 n 2 n i,j=1 κ(x i , x j ) -κ ⊤ K -1 κ. (8)

Constrained sparse formulation of the center

While the box constraint on the coefficients requires advanced optimization techniques, it is easy to satisfy the equality constraint (see ( 2)). The constrained optimization problem becomes α eq = arg min α1...,αn

1 n n i=1 Φ(x i ) - i∈I α i Φ(x i ) 2 subject to 1 ⊤ α eq = 1,
where 1 is a column vector of 1's. By using the Lagrangian multipliers, we obtain

α eq = α - K -1 1(1 ⊤ α -1) 1 ⊤ K -1 1 , ( 9
)
where α is the unconstrained solution, as given in [START_REF] Kemmler | One-class classification with gaussian processes[END_REF]. One may also include a regularization term, as above.

SOME THEORETICAL RESULTS

Independently of the algorithm, one is considering a set of samples in order to estimate the true expectation. Let

c ∞ = X Φ(x) dP (x)
be the true expectation, where P (x) is the probability distribution generating the samples x 1 , x 2 , . . . , x n . From these samples, one can give an estimate of c ∞ by using the empirical first moment c n , as defined in expression (4). The accuracy of such approximation is

ǫ 0 = c n -c ∞ .
Based on the Hoeffding's inequality, it is shown in [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF] (see also [START_REF] Bousquet | Introduction to Statistical Learning Theory[END_REF]) that, with probability at least 1 -δ over the choice of a random set of n samples, we have

nǫ 2 0 ≤ sup x∈X κ(x, x) 2 + √ -2 ln δ 2 .
By the symmetry of the i.i.d assumption, we can bound the probability that a new sample x, generated from the same probability distribution, is beyond the boundary defined by a one-class classification method, as given by the following proposition: Proposition 1. Consider the sphere centered on c I with radius max i=1,...,n Φ(x i )c I + 2ǫ 0 + 2 c nc I . Then, with probability at least 1 -δ over the choice of a random set of n samples, we can bound the probability that a new sample x is outside this sphere, with

P Φ(x) -c I > max i=1,...,n Φ(x i ) -c I + 2ǫ 0 + 2 c n -c I ≤ 1 n + 1 .
Proof. To show this, we consider Φ(x)c I and apply the triangle inequality twice, we get

Φ(x) -c I ≤ Φ(x) -c n + c n -c I ≤ Φ(x) -c ∞ + ǫ 0 + c n -c I ,
where the first inequality follows from approximating the fullmodel center by a subset of samples, while the second inequality from estimating the expected center by a finite set of n samples. Equivalently, we have for any x i :

Φ(x i ) -c I ≥ Φ(x i ) -c n -c n -c I ≥ Φ(x i ) -c ∞ -ǫ 0 -c n -c I .
Therefore, we get

P Φ(x) -c I > max i=1,...,n Φ(x i ) -c I + 2ǫ 0 + 2 c n -c I ≤ P Φ(x) -c ∞ > max i=1,...,n Φ(x i ) -c ∞ ≤ 1 n + 1
where the first inequality follows from the above inequalities, and the last inequality is due to the symmetry of the i.i.d assumption, considering n + 1 samples drawn from the same distribution.

As a special case of this proposition, consider the fullmodel for the empirical center, namely I = {1, 2, . . . , n}. In this case we get the relation given in [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF]Chapter 5]:

P Φ(x) -c n > max i=1,...,n Φ(x i ) -c n + 2ǫ 0 ≤ 1 n + 1 .
We can extend this result to the solution defined by considering that the samples defined by indices I are outliers, thus not inside the sphere. The following proposition can be easily proven using the same steps as in the proof of Proposition 1.

Proposition 2. Consider the same setting as in Proposition 1, where the indices in I define the outliers with |I| the number of outliers. Then, with probability at least 1 -δ over the choice of a random set of n samples, we can bound the probability that a new sample x is outside the sphere excluding outliers, with

P Φ(x) -c I > min i∈I Φ(x i ) -c I + 2ǫ 0 + 2 c n -c I ≤ |I| n + 1 .
It is worth noting that in both propositions, the error c nc I is minimized by our approach, as given by (8).

EXPERIMENTS

Since there are few benchmark datasets for one-class classification methods, multiclass tasks are often considered. A multiclass classification task can be tackled by using oneclass machines: each class is defined by a one-class classifier, and subsequently we get the decision rule by combining these classifiers. In practice, the model parameters are estimated by considering a subset of the target class (n train samples), and tested over the remaining samples (n test ), some from the target class and all the samples from the other classes.

To illustrate the relevance of the proposed approach, we have tested the proposed methods on real datasets well-known in the literature of one-class machines [START_REF] Wang | Structured one class classification[END_REF]: the IRIS dataset with 150 samples in 3 classes and 4 features, and WINE with 178 samples in 3 classes and 13 features. These datasets are available from the UCI machine learning repository. For experiments on IRIS data, we used only third and fourth features, as often investigated.

The Gaussian kernel was applied, with κ(x i , x j ) = exp( x ix j 2 /2σ 2 ). To estimate the classification error, a ten-fold cross-validation was used, with parameters optimized by a grid search over ν ∈ {2 -5 ; 2 -4 ; • • • ; 2 4 ; 2 5 } and σ ∈ {2 -5 ; 2 -4 ; • • • ; 2 4 ; 2 5 }. In order to provide comparable results, the number of SVs was fixed for all methods, by considering the optimal configuration for the classical one-class SVM. Table [START_REF] Zhang | Adaptive and online one-class support vector machine-based outlier detection techniques for wireless sensor networks[END_REF] gives the results in terms of classification error for each of the proposed methods, and compared to the classical one-class SVM. We also included the ratio of common SVs between the latter and each of the proposed methods, as well as the mean values. The computational cost of these machines, using the best configuration, are illustrated in terms of CPU time, as estimated on a 64-bit Matlab running on a Mackbook Pro with a 2.53 GHz Intel Core 2 Duo processor and 4 GB RAM.

CONCLUSION

In this paper, we studied the problem of one-class classification. By offering three one-class classification methods, we have shown that we can achieve a classification one-class while minimizing the classification error and especially with less computing time. The relevance of our approach is illustrated by experiments on well-known datasets. In future works, we study an online strategy for one-class classification, as well as other sparsification rules such as the coherence criterion. 1. Experimental results with the classification error for each one-class classifier, and the total 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.

Classical