Chafic Saidé

Régis Lengellé

Paul Honeine
email: paul.honeine@utt.fr

Cédric Richard

Roger Achkar

Chafic Saide

Régis Lengelle

Nonlinear Adaptive Filtering using Kernel-based Algorithms with Dictionary Adaptation

Keywords: Adaptive filters, machine learning, nonlinear systems, online system identification, kernel methods

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Nonlinear system modeling and identification have been widely studied in the last decades [START_REF] Aihara | Adaptive filtering for stochastic risk premia in bond market[END_REF][START_REF] Yu | Adaptive filtering with adaptive p-power error criterion[END_REF], using for example Volterra filters [START_REF] Schetzen | The Volterra and Wiener theories of nonlinear systems[END_REF][START_REF] Wiener | Nonlinear Problems in Random Theory[END_REF][START_REF] Ogunfunmi | Adaptive Nonlinear System Identification: The Volterra and Wiener Model Approaches[END_REF] and neural networks [START_REF] Haykin | Neural Networks and Learning Machines[END_REF]. More recently, function approximation based on Reproducing Kernel Hilbert Spaces (RKHS) defined 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][START_REF] Honeine | On-line nonlinear sparse approximation of functions[END_REF]. Kernel-based regression generally extends standard methods used in the linear case by replacing the computation of some inner products by a Mercer kernel [START_REF] Mercer | Functions of positive and negative type and their connection with the theory of integral equations[END_REF]. This provides an efficient way to implicitly map the data into a high, even infinite, dimensional RKHS and apply the original algorithm in this space without making direct reference to the nonlinear mapping of input vectors [START_REF] Aronszajn | Theory of reproducing kernels[END_REF][START_REF] Aizerman | Theoretical foundations of the potential function method in pattern recognition learning[END_REF]. Kernel-based regression requires the use of the Gram matrix composed of all the inner products (evaluated in the RKHS) of the different input vectors. Accordingly, the size of the Gram matrix linearly increases with time. Without introducing some rules to control the model complexity, this makes these methods intractable for real time applications.

Various algorithms were proposed to overcome this difficulty. The basic idea consists in introducing a new input sample to the model if it contributes in reducing the resulting approximation error and of removing the less relevant element, if necessary. A sparsification rule based on orthogonal projection was proposed in [START_REF] Dodd | A theory of iterative sparse solutions to strict interpolation in reproducing kernel hilbert spaces[END_REF][START_REF] Dodd | Function estimation in hilbert space using sequential projections[END_REF][START_REF] Dodd | Sparse,online learning in reproducing kernel hilbert spaces[END_REF][START_REF] Phonphitakchai | Stochastic meta descent in online kernel methods[END_REF]. The computational burden of this method is its main drawback. Another method considering the linear dependence criterion is detailed in [START_REF] Engel | Sparse online greedy support vector regression[END_REF][START_REF] Engel | Kernel recursive least squares[END_REF].

The coherence criterion, introduced in [START_REF] Honeine | Modélisation parcimonieuse non linéaire en ligne par une méthode á noyau reproduisant et un critère de cohérence[END_REF][START_REF] Richard | Online prediction of time series data with kernels[END_REF], has been used to control the number of kernel functions introduced in the dictionary † so as to limit the model complexity and the size of the Gram matrix. When using the coherence criterion, inclusion of a new input sample into the dictionary is performed if this element cannot be approximated by a linear combination (in the RKHS) of the actual elements of the dictionary. The main result in [START_REF] Richard | Online prediction of time series data with kernels[END_REF] is that the use of the coherence criterion induces a finite size of the dictionary when time goes to infinity.

However, for all previous approaches, each element introduced into the dictionary remains permanently unchanged, even if the non stationarity makes it later useless for estimating the solution. Therefore, in this paper we also consider the dictionary elements as adjustable model parameters in order to minimize the instantaneous output quadratic error. The main contribution of this paper is to perform dictionary adaptation. We propose an adaptation scheme based on a stochastic gradient algorithm. Our algorithm fulfills the constraints resulting from the coherence criterion. We also show that there is no universal dictionary adaptation algorithm, in the sense that each kernel function requires a specific study. We consider here the widely used gaussian, exponential and polynomial kernels and, depending on the used kernel function, we propose specific heuristics for selecting, at each iteration, the dictionary element(s) to be adapted.

In order to demonstrate the effectiveness of our method, we compare the performance obtained with and without dictionary adaptation on the following adaptive filtering algorithms: Kernel Affine Projection Algorithm (KAPA), Kernel Normalized Least Mean Squares (KNLMS), and Kernel Recursive Least Squares (KRLS). We analyze the performance obtained in terms of reduction of the approximation error and/or reduction of the size of the dictionary on several simulated and real data sets.

The paper is organized as follows. In section II, we propose a brief review of nonlinear system modeling methods that make the use of complexity control techniques. Investigations to adapt the elements of a dictionary using polynomial kernel functions and unit-norm kernel functions with radial basis functions are explored in section III. Section IV gives a brief review of the kernel-based adaptive algorithms used in our experiments (i.e., KAPA, KNLMS, and KRLS). The results obtained on synthetic and real data sets are presented in section V. Finally, our conclusion is presented in section VI.

A PRIMER ON KERNEL METHODS AND MODEL REDUCTION

Consider a nonlinear online prediction problem and let u n ∈ U ⊂ R ℓ be the input data vector and d n ∈ R the corresponding desired output of the model at time step n. Mapping the input data u n into a high-dimensional space using some nonlinear function ϕ(•) allows the use of linear modeling techniques to the transformed data ϕ(u n ).

Reproducing kernel Hilbert space and nonlinear regression

Considering a Hilbert space of real-valued functions ψ(•) denoted H, a compact U ⊂ R ℓ , and let

•, • H be the inner product in H. According to Riesz's representation theorem in [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], there exists a unique positive definite kernel denoted by κ(•, u j ) called representer of evaluation at u j , which satisfies: [START_REF] Mercer | Functions of positive and negative type and their connection with the theory of integral equations[END_REF], there exists a function ϕ(•) so that κ(u i , u j ) = ϕ(u i ), ϕ(u j ) H . The kernel then evaluates the inner product of any pair of elements of U mapped into H without the explicit knowledge of ϕ(•) or H. This is known as the kernel trick.

ψ(u j ) = ψ(•), κ(•, u j ) H , ∀ψ ∈ H, for every fixed u j ∈ U. Replacing ψ(•) by κ(•, u i ), yields κ(u i , u j ) = κ(•, u i ), κ(•, u j ) H , for all u i , u j ∈ U. If κ(•, •) is a Mercer kernel
Let κ(•, •) : U × U → R be a kernel and H the RKHS associated with κ(•, •). Using the kernel trick, linear algorithms expressed in term of inner products are easily transformed into nonlinear ones. The least-squares algorithm consists of determining a function ψ * (•) of H that minimizes a cost function i.e., the sum of the squared error between the actual response of the system (desired output) d j and the corresponding output of the model ψ(u j ) = ψ(•), κ(•, u j ) H , i.e., ψ * (•) is the solution of the regularized risk functional

min ψ∈H n j=1 |d j -ψ(u j )| 2 + ξ g( ψ 2 ), ( 1
)
where • is the norm in the RKHS, ξ g( ψ 2 ) is a regularization term (ξ ≥ 0) and g(•) is a strictly monotonic function on R + . According to the representer theorem [START_REF] Schölkopf | A generalized representer theorem[END_REF][START_REF] Kimeldorf | Some results on Tchebycheffian spline functions[END_REF][START_REF] Wahba | Spline models for observational data[END_REF], the optimum model ψ * (•) can be written

ψ * (•) ≡ ψ n (•) = n j=1 α j κ(•, u j ). (2)
The minimization of the risk functional (1) leads to an optimization problem with n variables which are the coefficients of the model {α j } j=1•••n .

A model reduction method using the coherence criterion

The main problem with online identification methods is the increasing number of observations with time so that adaptation of the model (2) for every new observation becomes more and more complex.

The sparsification is the process of selecting the most relevant kernel functions among the past observations and use them to predict the model. The selected kernel functions form the dictionary and the number of these functions denotes its order m. Many sparsification measures have been used to characterize a dictionary in linear sparse approximation techniques [START_REF] Engel | Kernel recursive least squares[END_REF][START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF][START_REF] Sun | On the size of the online kernel sparsification dictionary[END_REF][START_REF] Csat | Sparse representation for gaussian process models[END_REF]. It was also introduced by Mallat and Zhang for matching pursuit [START_REF] Mallat | Matching pursuit with time-frequency dictionaries[END_REF]. In [START_REF] Richard | Online prediction of time series data with kernels[END_REF], Richard, Bermudez and Honeine proposed to define the coherence parameter as

µ = max i =j | κ(•, u i ), κ(•, u j ) H | = max i =j |κ(u i , u j )|,
where κ is a unit-norm kernel. If κ is not a unit-norm kernel, κ(u i , u j ) can be replaced by κ(u i , u j )/ κ(u i , u i )κ(u j , u j ). The parameter µ is the largest absolute value of the off-diagonal entries in the Gram matrix and it reflects the largest cross-correlations in the dictionary. Therefore, coherence equals zero for every orthogonal basis.

This measure of linear dependence in H is used to decide if, at each time step n, a candidate kernel function κ(•, u n ) must be introduced into the dictionary or not. If this candidate kernel function can be represented by a linear combination of the kernel functions of the dictionary it is dropped, else it is included into the dictionary and the order m is incremented (m = m + 1). Introduction of the candidate function κ(•, u n ) into the dictionary is decided when coherence remains smaller than a threshold µ 0 ,

max wj ∈Jn-1 |κ(u n , u wj )| ≤ µ 0 , (3)
where µ 0 is a parameter ∈ [0, 1[ determining the level of sparsity and the coherence of the dictionary,

J n-1 denotes a subset of size m of {1, • • • , n -1} and D n-1 = {κ(•, u wj )} m j=1
is the dictionary at time step n -1. A dictionary fulfilling (3) is said µ 0 -coherent ‡ . At time step n, if κ(•, u n ) satisfies the condition (3) it will be introduced into the dictionary which becomes

D n = D n-1 ∪ {κ(•, u n )}
The key point in the use of the coherence criterion is, as shown in [START_REF] Richard | Online prediction of time series data with kernels[END_REF], that the size of the dictionary remains finite with time n so the obtained model becomes of fixed-size. It can be written

ψ n (•) = m j=1 α n,j κ(•, u wj ), (4)
where m ≡ m(n) converges to some finite value as n → ∞ for any value of µ 0 < 1 and

α n = [α n,1 , • • • , α n,m
] t is the optimal solution vector at time step n.

DICTIONARY ADAPTATION

In all previous approaches using sparsification criterion to select dictionary elements, a kernel function κ(•, u wi ) once introduced into the dictionary, remains unchanged even if it becomes useless in evaluating the model output. This is due to possible non-stationarity. In our method, the dictionary is adapted at each time step in order to minimize the instantaneous quadratic output error. In our approach, the {u wj } m j=1 become variables to be optimized, as well as {α n,j } m j=1 . The main constraint is that the dictionary must remain µ 0 -coherent after adaptation.

Let D n = {κ(•, u wi )} m
i=1 be a µ 0 -coherent dictionary at time step n and let m denotes its cardinality. Our objective is to adjust the elements of the dictionary D n to obtain D A n using an adaptive method. We use an alternate directions optimization scheme thus, in a first step, the optimal value of the solution vector

α n = [α n,1 , • • • , α n,
m ] t is found using an online adaptive algorithm. ‡ A dictionary is µ 0 -coherent if all its kernel functions verify the following condition: max i =j |κ(uw i , uw j )| ≤ µ 0 . Second, the dictionary is adapted with the objective of reducing the instantaneous quadratic output error e 2 n and in such a way that the new dictionary still fulfills the coherence criterion, as follows:

min uw i ∈Jn e 2 n , (5)
subject to the constraint

max i =j,i,j∈Jn |κ(u wi , u wj )| ≤ µ 0 ,
where

e n = d n -ψ n (u n ) = d n - m i=1 α n,i κ(u n , u wi ).
Here, d n represents the desired output of the model.

The global adaptation algorithm is represented in Figure 1. The i-th element of the dictionary u wi is updated according to:

u A wi = u wi -ν n g wi ∀i = 1 • • • m, (6)
where u A wi is the i-th dictionary element after adaptation, g wi is the gradient of the instantaneous quadratic error with respect to u wi , and ν n represents the gradient step size used to adjust the elements of the dictionary. Obviously, the gradient step size must be selected so as the new dictionary remains in the feasible domain.

• By construction of our algorithm, the coherence criterion is satisfied at each time step n before adaptation

max i =j |κ(u wi , u wj )| ≤ µ 0 , ∀i, j = 1 • • • m, (7)
and this must remain true after adaptation, so we must have

max i =j |κ(u A wi , u A wj )| ≤ µ 0 . ∀i, j = 1 • • • m. (8)
• Using ( 5), the gradient with respect to u wi is:

g wi = ∇ u w i (e 2 n ) = 2 e n ∇ u w i d n - m i=1 α n,i κ(u n , u wi ) = -2 e n α n,i ∇ u w i κ(u n , u wi ) . (9)
From the expression of the gradient g wi [START_REF] Honeine | On-line nonlinear sparse approximation of functions[END_REF], it is obvious that the choice of the kernel function to be used has a major influence in the adaptation process. In this paper, we explore the radial basis kernel, which belongs to the class of translation-invariant kernels, and the polynomial kernel.

Case of the radial basis kernel

Most of translation-invariant kernels, mainly the gaussian and the exponential kernel functions, can be expressed in the form:

κ(u i , u j ) = f ( u i -u j 2 ), (10)
where f ∈ C ∞ . A sufficient condition for this function to be a valid positive-definite kernel is the complete monotonicity of the function f , i.e.,

(-1) k f (k) (r) ≥ 0, ∀r ≥ 0, (11)
where

f (k) (•) denotes the k-th derivative of f (•) [29]
. From [START_REF] Mercer | Functions of positive and negative type and their connection with the theory of integral equations[END_REF], we get

∇ u j κ(u i , u j ) = -2(u i -u j ) f (1) ( u i -u j 2 ). ( 12
)
Note that, in the case of a radial basis kernel, the condition in [START_REF] Smola | A tutorial on support vector regression[END_REF] for any two elements of a µ 0coherent dictionary is expressed as follows:

f ( u wi -u wj 2 ) ≤ µ 0 . ( 13
)

Gradient and coherence condition for the gaussian kernel

The kernel function for the i-th element of the dictionary is:

κ(•, u wi ) = exp - • -u wi 2 2σ 2 ,
where σ is the bandwidth of the kernel. Taking the derivative of this kernel w.r.t. u wi , the gradient of the instantaneous quadratic error defined in ( 9) and ( 12) yields

g wi = -2 e n α n,i σ 2 κ(u n , u wi ) (u n -u wi ).
the coherence condition for the gaussian kernel, according to expressions [START_REF] Smola | A tutorial on support vector regression[END_REF] and ( 8) is:

u A wi -u A wj 2 ≥ -2 σ 2 ln(µ 0 ), with -2 σ 2 ln(µ 0 ) ≥ 0 because µ 0 ∈ [0, 1[.

Gradient and coherence condition for an exponential kernel

The kernel function for the i-th element of the dictionary is:

κ(•, u wi ) = exp - • -u wi 2σ 2 ,
where σ is the bandwidth of the kernel function. The gradient of the instantaneous quadratic error with respect to u wi at time step n is:

g wi = - e n α n,i σ 2 u n -u wi κ(u n , u wi ) (u n -u wi ).
Obviously, from the above expression, u nu wi must be non zero, meaning that when introducing u n in the dictionary, the newly introduced element cannot be adjusted.

The coherence condition for an exponential kernel, according to equations ( 7) and ( 8) is:

u A wi -u A wj 2 ≥ -2σ 2 ln(µ 0 ) 2 .

Finding a feasible gradient step size in the case of a radial basis kernel

The main issue is clearly to find an adequate gradient step size to simultaneously adjust all the dictionary elements without violating the coherence criterion. For any pair of dictionary elements, we obtain

u A wi -u A wj = δu -ν n δg ∀i, j = 1 • • • m,
where δu = u wi -u wj and δg = g wi -g wj . Because of the constraint induced by the coherence criterion, ν n cannot be chosen arbitrarily. The problem is now to determine an appropriate ν n at time step n in order to adapt the dictionary. The coherence condition (8) leads to

f ( δu -ν n δg 2 ) ≤ µ 0 . (14)
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 (14) becomes

-2 δg 2 ν 2 n -2ν n δu t δg)f (1) ( δu 2 + µ 0 -f ( δu 2 ) ≥ 0 -2 δg 2 f (1) δu 2 ν 2 n + 2δu t δgf (1) ( δu 2 ν n + µ 0 -f ( δu 2 ) ≥ 0 (15)
It is very important to note that the first derivative of a radial basis function is always negative on R + , which means that the coefficient of the second degree term in ν n is always positive, namely

-2 δg 2 f (1) ( δu 2 ) > 0.
Second, from [START_REF] Dodd | A theory of iterative sparse solutions to strict interpolation in reproducing kernel hilbert spaces[END_REF] we have µ 0 -f ( δu 2 ) ≥ 0 as the coherence criterion is fulfilled at any time step.

The reduced discriminant of this quadratic expression in ν n is:

∆ = δu t δg f (1) ( δu 2 ) 2 + 2 δg 2 f (1) ( δu 2 )(µ 0 -f ( δu 2 ).
Let ν ij n be a valid gradient step size for adapting the two dictionary elements u wi and u wj (at time step n). If ∆ < 0, there is no constraint on the value of the step size ν ij n . Otherwise, if ∆ ≥ 0, the boundary points of (15) will be ν i,j-and ν i,j+ as follows:

ν i,j± = -δu t δg f (1) ( δu 2 ) ± √ ∆ -δg 2 f (1) ( δu 2 ) , 0 0 0 0 ν i,j- ν i,j+ ν 0 ν ij n = ν 0 ν i,j- ν i,j+ ν 0 ν ij n = ν 0 ν i,j- ν i,j+ ν 0 ν ij n = ν i,j- ν i,j- ν i,j+ ν 0 ν ij n = ν i,j- Figure 2.
An illustration showing all the possibilities to select the gradient step size between any two dictionary elements uw i and uw j in the case of a radial basis kernel. ν 0 > 0 is a reference step size. and the interval of possible values for ν ij n will be ] -∞, ν i,j-] ∪ [ν i,j+ , +∞[, since the quadratic expression in (15) must be positive:

ν i,j- ν i,j+ + - +
Figure 2 shows all the possibilities to select the gradient step size when adapting any two dictionary elements u wi and u wj in the case of a radial basis kernel. Obviously, for each pair of dictionary elements (u wi , u wj ), ν ij n = 0 always belong to the feasible domain, since there is no adaptation and the initial dictionary was µ 0 -coherent. This is also the reason for which we can always use Taylor series around ν n = 0. Finally because the coefficient of the second degree term in ν n and the constant term in [START_REF] Dodd | Sparse,online learning in reproducing kernel hilbert spaces[END_REF] are positive, ν i,j-and ν i,j+ have the same sign.

Interpretation of the two bounds for ν n is straightforward. When adjusting any two dictionary elements, each one according to its own gradient of the quadratic error, the interval ]ν i,j-, ν i,j+ [ must be eliminated to avoid any overlap of the influence regions of the two considered elements, and thus violation of the coherence constraint [START_REF] Suykens | Least squares support vector machines[END_REF]; see Figure 3. We initially select a reference step size ν 0 > 0, similarly as for all adaptive algorithms with a fixed step size. After calculating the m(m -1)/2 intervals [ν i,j-, ν i,j+ ], the single adaptation step ν n used for all dictionary elements is selected as illustrated in Figure 2:

9 ν i,j+ ν i,j- u A wi u A wi u wi u A wj -g wi
• if max i,j ν i,j+ ≤ 0 ⇒ ν 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 ; • 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 the positive part of ν i,j-. 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. Other heuristics for dictionary adaptation could be considered, but our main objective is only to illustrate the potential efficiency of dictionary adaptation.

Case of a Polynomial Kernel

Gradient and coherence criterion

The polynomial kernel function is of the form:

k(u i , u j ) = f (u t i u j ) = (1 + u t i u j ) β , (16)
with β > 0. This kernel function is not unit-norm, hence, because of the use of the coherence criterion, we need to normalize it to become a unit-norm kernel [START_REF] Richard | Online prediction of time series data with kernels[END_REF]. This is achieved by substituting k(u i , u j ) with:

k(u i , u j ) k(u i , u i ) k(u j , u j ) .
Therefore, equation [START_REF] Ogunfunmi | Adaptive Nonlinear System Identification: The Volterra and Wiener Model Approaches[END_REF] becomes:

e n = d n - m i=1 α n,i k(u n , u wi ) k(u n , u n ) k(u wi , u wi ) . ( 17
)
The gradient g wi of the instantaneous quadratic error with respect to u wi is:

g wi = -2 β e n α n,i k(u n , u wi ) k(u n , u n ) k(u wi , u wi ) u n k * (u n , u wi ) - 1 2 u wi k * (u wi , u wi ) , ( 18
)
where k * (u i , u j ) = (1 + u t i u j ). The coherence criterion [START_REF] Smola | A tutorial on support vector regression[END_REF], when applied to the polynomial kernel, yields:

|(1 + u t wi u wj ) β | (1 + u t wi u wi ) β (1 + u t wj u wj ) β ≤ µ 0 , (19)
which can be written:

(1 + u t wi u wj ) 2 (1 + u wi 2 ) (1 + u wj 2 ) ≤ (µ 0 ) 2/β . (20)
Considering [START_REF] Richard | Online prediction of time series data with kernels[END_REF], it seems to be challenging to adjust simultaneously all the dictionary elements while fulfilling the coherence criterion (3). This explains our choice of selecting, for the polynomial kernel, only one element of the dictionary to be adapted at each time step. Heuristics for selecting this element are presented in section 3.2.4.

Search for possible values of the gradient step size

We must now be sure that, after adaptation of the selected element, the new dictionary still fulfills the coherence criterion. Let u wi denote the dictionary element to be adapted. Introducing [START_REF] Haykin | Neural Networks and Learning Machines[END_REF] in [START_REF] Richard | Online prediction of time series data with kernels[END_REF] gives:

1 + (u A wi ) t u wj 2 (1 + u A wi 2 ) (1 + u wj 2 ) ≤ (µ 0 ) 2/β ∀j = i, j = 1 • • • m. (21)
Since (u A wi ) t u wj = (u wi -ν n g wi ) t u wj = u t wi u wj -ν n u t wj g wi , and

u A wi 2 = u wi 2 -2ν n u t wi g wi + ν 2 n g wi 2 ,
inequality (21) becomes:

k * 2 ij -2 ν n u t wj g wi k * ij + ν 2 n u t wj g wi 2 k * jj k * ii -2ν n u t wi g wi + ν 2 n g wi 2 ≤ (µ 0 ) 2/β , (22) 11
where k * ij = k * (u wi , u wj ) = (1 + u t wi u wj ). Setting D = k * jj (µ 0 ) 2/β , ( 22) becomes:

u t wj g wi 2 -D g wi 2 ν 2 n -2 k * ij u t wj g wi -Du t wi g wi ν n + k * 2 ij -Dk * ii ≤ 0, (23)
which is of the form:

A ν 2 n + 2B ν n + C ≥ 0, (24)
with

A = D g wi 2 -u t wj g wi 2 , B = k * ij u t wj g wi -Du t wi g wi , C = Dk * ii -k * 2 ij .
Here, C is always positive for a coherent dictionary, this results from the fulfillment of the coherence criterion:

k ij k ii k jj ≤ µ 0 ⇒ k * 2 ij k * ii k * jj ≤ (µ 0 ) 2/β .
The reduced discriminant is:

∆ = B 2 -AC = k * ij u t wj g wi -Du t wi g wi 2 -D g wi 2 -(u t wj g wi ) 2 Dk * ii -k * 2 ij .
As can be seen from ( 24), selecting a valid value of ν n depends on the sign of A. The following analysis, based on the study of the positiveness of a second degree polynomial, considers all possible cases.

• Case of A < 0. In this case ∆ ≥ 0 if ∆ > 0 then the polynomial defined in [START_REF] Kimeldorf | Some results on Tchebycheffian spline functions[END_REF] has two roots with opposite signs:

ν i,j ± = -k * ij u t wj g wi -Du t wi g wi ± √ ∆ D g wi 2 -u t wj g wi 2 ,
with:

ν i,j- ν i,j+ - + - ν n ∈ [ν i,j-, ν i,j+ ] with ν i,j-≤ 0 ≤ ν i,j+
, because ν n = 0 is always possible since the dictionary is not adapted in this case.

if ∆ = 0, then B 2 = AC and since A < 0 and C ≥ 0, this implies that B = 0 and ν i,j-= ν i,j+ = -B A = 0.

• Case of A > 0. In this case, three possibilities must be considered:

if ∆ > 0, the two roots have the same sign since C A ≥ 0, with: if ∆ < 0 the polynomial does not have roots but is always positive, meaning that there is no constraint for selecting ν n ∈ ] -∞, +∞[.

ν i,j- ν i,j+ + - + and ν n ∈ ] -∞, ν i,j-] ∪ [ν i,
if ∆ = 0, the polynomial is positive or null and there is no constraint for selecting ν n .

• Case of A = 0. The unique solution is

ν i,j = -C 2B -if B > 0, this means that ν n ∈ ]ν i,j , +∞ [. -if B ≤ 0, then -C 2B ≥ 0 and ν n ∈ ] -∞, ν i,j [.

Selection of the gradient step size

After selecting the element u wi to adapt, the calculation is made to satisfy the (m -1) constraints with the other dictionary elements. We first select, as usual, a reference step size ν 0 > 0. Second, for j = i and j = 1 • • • m, the value of ν ijn corresponding to the couple (u wi , u wj ) is selected as follows:

1. Case of A > 0 and ∆ > 0

• if ν i,j-≤ ν i,j+ ≤ 0 ⇒ ν ijn = ν 0 • if 0 ≤ ν i,j-≤ ν 0 ≤ ν i,j+ ⇒ ν ijn = ν i,j- • if 0 ≤ ν i,j-≤ ν i,j+ ≤ ν 0 ⇒ ν ijn = ν 0 or ν ijn = ν i,j- • if 0 ≤ ν 0 ≤ ν i,j-≤ ν i,j+ ⇒ ν ijn = ν 0 or ν ijn = ν i,j- 2. Case of A > 0 and ∆ ≤ 0 • if such is the case, we select ν ijn = ν 0 3. Case of A < 0 and ∆ > 0 in this case ν i,j-≤ 0 ≤ ν i,j+ • if ν i,j+ ≤ ν 0 ⇒ ν ijn = ν i,j+ • if ν 0 ≤ ν i,j+ ⇒ ν ijn = ν 0 or ν ijn = ν i,j+
4. Case of A < 0 and ∆ = 0

• we have ν ijn = 0 13 5. Case of A = 0 • if -C 2B < 0 ⇒ ν ijn = ν 0 . • if 0 ≤ -C 2B ≤ ν 0 ⇒ ν ijn = -C 2B . • if ν 0 ≤ -C 2B ⇒ ν ijn = ν 0 .
When the (m -1) values ν ijn are obtained, we select ν n as ν n = min j =i, j=1•••m ν ijn .

Selection of the element to be adapted

Several approaches can be adopted to select the element that will be adapted. We only investigated two methods. The first one chooses the element that corresponds to the smallest absolute value of the components of the solution vector

α n = [α n,1 , • • • , α n,m ] t .
The second consists of calculating the partial gradient of the quadratic error with respect to each element of the dictionary. Then we select the one with the highest norm of the gradient. The advantage of the first heuristic is its simplicity, while the second gives better results but requires more computations. Both methods were implemented and are compared in the experimentation section.

A BRIEF REVIEW OF SOME ONLINE LEARNING ALGORITHMS

In this paper, and for reasons only related to the possibility to compare the performance obtained with previously published results, we use the Kernel Recursive Least Squares algorithm (KRLS), the Kernel Affine Projection Algorithm (KAPA), and the Kernel Normalized Least Squares algorithm (KNLMS). These algorithms are used with and without dictionary adaptation. The output of the model is:

ψ n (•) = m j=1 α n,j κ(•, u wj ). (25)
The optimal solution vector α n is found using one of the above mentioned online adaptive algorithm.

Kernel Recursive Least Squares algorithm (KRLS)

In the Kernel Recursive Least Squares algorithm, the optimization problem (1) can be written as follows [START_REF] Honeine | Modélisation parcimonieuse non linéaire en ligne par une méthode á noyau reproduisant et un critère de cohérence[END_REF]:

min α d n -H n α 2 + ξα t K n α, (26)
where K n is the m × m Gram matrix of the dictionary at time n, ξ is a regularization coefficient, and H n is an (n × m) matrix with (i, j)-th element is κ(u i , u wj ). At time step n the solution of the problem (26) will be [START_REF] Honeine | Modélisation parcimonieuse non linéaire en ligne par une méthode á noyau reproduisant et un critère de cohérence[END_REF]:

α n = P n H t n d n , (27)
provided that P n = (H n H t n + ξK n ) -1 exists. At time step n+1, there is a new input u n+1 to the model and one of the following cases can occur:

1. First case: max j=1•••m |κ(u n+1 , u wj )| > µ 0
The new input can be represented by the actual elements of the dictionary. It is not introduced in the dictionary. The matrix H n is updated by adding a (n + 1) -th row

h n+1 = [κ(u n+1 , u w1 ), • • • , κ(u n+1 , u wm )]
t and the new desired output d n+1 is added to the output vector d n . Update of the model is performed as follows:

α n+1 = α n + P n h n+1 1 + h t n+1 P n h n+1 d n+1 -h t n+1 α n , (28)
and

P n+1 = P n - P n h n+1 h t n+1 P n 1 + h t n+1 P n h n+1 . ( 29
) 2. Second case: max j=1•••m |κ(u n+1 , u wj )| ≤ µ 0
The new input cannot be represented by any actual element of the dictionary. Thus, it must be introduced in the dictionary. The order of the dictionary becomes (m + 1). First, the vector α n and P n are updated according to ( 28) and ( 29), leading to αn+1 and Pn+1 respectively.

Second, in order to introduce a new coefficient in the model, the new α n+1 and P n+1 are given by:

α n+1 =    αn+1 0    + d n+1 -h t n+1 α n 1 -h t n+1 Pn+1 h n+1    Pn+1 h n+1 1/h 0    , (30)
P n+1 =    Pn+1 0 n 0 t n 0    + 1 1 -h t n+1 Pn+1 h n+1 ×    -Pn+1 h n+1 1/h 0    -( Pn+1 h n+1 ) t 1/h 0 , (31)
where h 0 = κ(u n+1 , u n+1 ).

Kernel Affine Projection Algorithm

The Kernel Affine Projection Algorithm, studied in the linear case in [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] Esfand Abadi | A family of set-membership affine projection adaptive filter algorithms[END_REF] and in the Kernelized APA in [START_REF] Liu | Kernel affine projection algorithms[END_REF][START_REF] Liu | Kernel Adaptive Filtering: A Comprehensive Introduction[END_REF], results from the following optimization problem (while neglecting the regularized term in (1)):

min α d n -H n α 2 . ( 32
)
The optimal solution α n of the problem ( 32) is:

α n = (H t n H n ) -1 H t n d n , (33)
provided that (H t n H n ) -1 exists and H n is a n-by-m matrix whose (i, j)-th entry is κ(u n-i+1 , u wj ), and the i-th element of the column vector d n is d n-i+1 . At each time step n, the affine problem will be:

min α α -α n-1 2 subject to d n = H n α, (34)
where, at each time step n, only the p most recent inputs {u n , • • • , u n-p+1 } and observations

{d n , • • • , d n-p+1
} are used, so the i-th element of the column vector

d n is d n-i+1 , i = 1 • • • , p,
and H n is the p × m matrix whose (i, j)-th entry is κ(u n-i+1 , u wj ). This means that α n is obtained by projecting α n-1 onto the intersection of the p manifolds.

At time step n + 1, a new input is fed to the model. One of the two following possibilities occurs:

1. First case: max j=1•••m |κ(u n+1 , u wj )| > µ 0
In this case, κ(•, u n+1 ) is not introduced in the dictionary D n and the solution vector α n+1 is updated according to [START_REF] Richard | Online prediction of time series data with kernels[END_REF]:

α n+1 = α n + η H t n+1 (ǫ I + H n+1 H t n+1 ) -1 (d n+1 -H n+1 α n ), ( 35
)
where η is called the step-size control parameter and ǫ I is a regularization factor.

Second case:

max j=1•••m |κ(u n+1 , u wj )| ≤ µ 0
In this case, κ(•, u n+1 ) cannot be represented by the actual kernel functions of the dictionary D n and it is introduced as a new element and denoted κ(•, u wm+1 ). The resolution of [START_REF] Liu | Kernel Adaptive Filtering: A Comprehensive Introduction[END_REF] leads to the following recursive solution:

α n+1 =    α n 0    + η H t n+1   ǫ I + H n+1 H t n+1 ) -1 (d n+1 -H n+1    α n 0       . ( 36
)
The set of recursions ( 35) and ( 36) are called the Kernel Affine Projection Algorithm (KAPA).

Kernel Normalized Least Mean Squares algorithm (KNLMS)

This algorithm is similar to the Kernel Affine Projection Algorithm with a number of manifolds p equal to 1. In this case, the constraint appearing in problem (34) becomes d n = h t n α n where h n is a column vector with i-th entry equal to κ(u n , u wi ). Relations ( 35) and (36) become:

1. First case: max j=1•••m |κ(u n+1 , u wj )| > µ 0 α n+1 = α n + η ǫ + ||h n+1 || 2 d n+1 -h t n+1 α n h n+1 . (37)
2. Second case:

max j=1•••m |κ(u n+1 , u wj )| ≤ µ 0 α n+1 =    α n 0    + η ǫ + ||h n+1 || 2   dn+1 -h t n+1    α n 0       h n+1 , (38)
with

h n+1 = [κ(u n+1 , u w1 ), • • • , κ(u n+1 , u wm+1 )] t .

EXPERIMENTATIONS

In this section, we illustrate the benefits of dictionary adaptation. To compare results, we choose ν 0 , the reference step-size for dictionary adaptation, and µ 0 , the coherence parameter, in such a way that the final mean quadratic errors are the same, with and without adaptation, and we compare the final sizes of the dictionaries. Another possibility is to select ν 0 and µ 0 to get almost identical sizes of the dictionaries so we can directly compare the residual quadratic error.

The Kernel Recursive Least Squares Algorithm (KRLS) performance

Our first experiment consists in predicting the Henon time series generated by:

d n = 1 -γ 1 d 2 n-1 + γ 2 d n-2 , with d 0 = -0.3, d 1 = 0, γ 1 = 1.4, and γ 2 = 0.3.
Modeling is performed with a model of the form

ψ n (d n-1 , d n-2 )
and makes the use of a series of 2000 samples. The gaussian kernel is used with parameter σ = 0.35. To compare our results with the ones in [START_REF] Honeine | Modélisation parcimonieuse non linéaire en ligne par une méthode á noyau reproduisant et un critère de cohérence[END_REF], we use the same parameters settings and a coherence criterion µ 0 = 0.6. The adaptation step size is ν 0 = 0.05. Modeling is performed with and without dictionary adaptation (see Figure 4 and Figure 5).

With adaptation, the dictionary size decreases from 17 to 16 elements which corresponds to a decrease of 5.88%. The instantaneous quadratic error is averaged over the last 500 samples.

Without adaptation, the mean of the quadratic error is 0.01036 versus 0.001628 with adaptation.

This corresponds to a decrease of 84.29% of this error.

Kernel Affine Projection Algorithm (KAPA)

We now present two other experiments using KAPA algorithm. The first consists in predicting the logistic time series using gaussian and exponential kernels. The second is extracted from [START_REF] Richard | Online prediction of time series data with kernels[END_REF] and makes the use of the polynomial kernel.

Radial basis kernels

In this experiment, the radial basis kernel is used with the Logistic Map which is a deterministic chaotic non linear model of the form: with initial conditions a = 4 and d 0 = 0.2027. A zero mean white gaussian noise with standard deviation 0.25 is added to the desired output. Here again the experiment is conducted with and without dictionary adaptation. Note that in the case of radial basis kernels we adapt all the dictionary element using the same gradient step size.

d n+1 = ad n (1 -d n ), -1.5 -1 -0.5 0 0.5 1 1.5 -1.5
The coherence parameter is µ 0 = 0.3. Results using the gaussian kernel with a bandwidth 0.418, are shown in Figure 6. These results demonstrate that dictionary adaptation (using ν 0 = 0.2) leads to a strong reduction of the mean squared error which is 0.0016778 with adaptation versus 0.0060538 without adaptation. This corresponds to a reduction of 72.285% of the error, with the same final dictionary size m = 2.

Results using an exponential kernel with a bandwidth 0.5, are shown in Figure 7. They indicate that dictionary adaptation (ν 0 = 0.01) leads to a mean squared error of 0.0099111 with a dictionary size m = 2. Without adaptation, the mean squared error is 0.01601 with the same dictionary size m = 2. This means that with the same dictionary size m = 2, the quadratic error is reduced with 38.09%.

Polynomial kernel

In this experiment,we use the same benchmark parameters setting as in [START_REF] Richard | Online prediction of time series data with kernels[END_REF]. It consists of a one The KAPA algorithm is used with the following parameters setting: number of manifolds p = 3, step-size control parameter η = 0.01, regularization factor ǫ = 0.07. A set of 200 time series of 3000 samples each is used to compare KAPA and AKAPA [START_REF] Saide | Dictionary adaptation for online prediction of time series data with kernels[END_REF] using the Normalized Mean Squared Error (NMSE) which is computed over the last 500 samples according to:

NMSE = E      3000 i=2501 (d n -ψ n (u n , u n-1 )) 2 3000 i=2501 d 2 n      .
The expected value is estimated over the 200 realizations of the time series.

We also compute the final size of the dictionary m averaged over the 200 realizations. The dictionary adaptation is performed for the element giving the highest value of the norm of the gradient along with a polynomial kernel of degree β = 2 (see Figure 8). Also, the other strategy to choose the element to adapt which is the one giving the smallest value of |α n,i | was test with a polynomial kernel of degree β = 3 (see Figure 9).

Using Table I, the following observations can be extracted: • for the same average final dictionary size, the NMSE decreases of 50.95% if the adjusted element is the one with the highest gradient norm (β = 2), and of 17.44% if the adjusted element is the one with the the smallest value of |α n,i | (β = 2).

• for approximately the same NMSE, the average dictionary size has decreased of 26.54% if the adjusted element is the one with the highest gradient norm (β = 2), and of 25.05% if the adjusted element is the one with the the smallest value of |α n,i | (β = 2).

Kernel Normalized Least Mean Squares algorithm (KNLMS)

In this section we consider the same nonlinear benchmark as proposed in [START_REF] Dodd | Function estimation in hilbert space using sequential projections[END_REF]. The desired system output d n is given by:

d n = (0.8 -0.5e -d 2 n-1 )d n-1 -(0.3 -0.9e -d 2 n-1 )d n-2 + 0.1 sin(d n-1 π).
The initial conditions are (0.1;0.1). The output is corrupted with an additive zero mean white noise with standard deviation equal to 0.1. The model is of the form ψ n (d n-1 , d n-2 ).

The same parameters as in [START_REF] Richard | Online prediction of time series data with kernels[END_REF] • When performing dictionary adaptation, and for the same average final size of the dictionary, NMSE has decreased of 75.148% for the gaussian kernel and of 73.885% for the exponential kernel.

• For the same final average NMSE, the average size of the dictionary has decreased of 34.50%

for the gaussian kernel and of 31.54% for the exponential kernel.

Note that for both kernels, algorithm converges faster when performing dictionary adaptation, which is not surprising. KNLMS and KRLS are compared in Figure 12.

Sun spots prediction using KAPA

We now apply our adaptation algorithm to a real time series which consists of the number of sunspots per month between January 1749 and February 2012. This data set has been obtained from [START_REF] Hathaway | The sunspot cycle[END_REF]. For this time series, we use the Kernel Affine Projection Algorithm coupled with the gaussian kernel function. The model is of the form ψ n (u n-1 , u n-2 , u n-3 ). Using the following parameter settings η = 0.1, ǫ = 0.07, p = 3, the NMSE is calculated over the last 300 samples.

Comparing NMSE while obtaining approximately the same final dictionary size m ≈ 536, without adaptation, setting µ 0 = 0.12415 the NMSE is 0.016839 while with adaptation and setting µ 0 = 0.1 and ν 0 = 0.0175, we observe that the NMSE has decreased to 0.0033112 inducing a reduction of 80.336%. The details are shown in Figure 13. However, even in a non stationary environment, previous work never considered the dictionary as a set of parameters that could be adapted in order to further reduce the instantaneous quadratic output error.

In this paper, we have introduced a method for optimizing the dictionary elements on which the model is decomposed, simultaneously with the standard model parameters adaptation. The coherence was the sparsification criterion used. Our gradient-based method allowed iterative fulfillment of the coherence criterion used both for introducing a new element in the dictionary and for adapting the existing dictionary elements. We have presented the complete updating equations in the case of unit-norm gaussian and exponential kernels and the non unit-norm polynomial kernel.

We have also shown that there is no universal dictionary updating algorithm, which means that every kernel type requires a specific study. We have compared the results obtained with and without dictionary adaptation, both on simulated and real world data, and we have shown that, using dictionary adaptation, the estimation error can be dramatically reduced (for a similar final size of the dictionary), or the final size of the dictionary can be reduced (for a similar residual estimation error). Obviously, any other compromise can be selected between these two extreme possibilities.

Finally, among the possible extensions of this work, we are working on reducing the computational complexity which is strongly related to the (time varying) number of dictionary elements and we are considering different sparsification criteria such as linear dependence measures.

Figure 1 .

 1 Figure 1. Global adaptation algorithm.

Figure 3 .

 3 Figure 3. A 2D illustration showing the constraint of choosing ν ij n ≤ ν i,j-or ν ij n ≥ ν i,j+ to avoid the overlap of the influence regions of u A wi and u A wj .

  j+ , +∞[. Of course ν n = 0 lies in the interval of possible values.

Figure 4 .

 4 Figure 4. Henon system modeling using KRLS with µ 0 = 0.6 and ν 0 = 0.05.

Figure 5 .µ 0

 50 Figure 5. Quadratic error of Henon system modeling using KRLS with µ 0 = 0.6 and ν 0 = 0.05.

Figure 6 .

 6 Figure 6. Logistic Map: Quadratic Error -KAPA with gaussian kernel

Figure 7 .d n = s 2 n

 72 Figure 7. Logistic Map: Quadratic Error -KAPA with exponential kernel

  The adapted element is the one with the highest gradient norm, β = 2 µ 0 = 0.3 m=1.62, NMSE=0.27649 µ 0 = 0.15, ν 0 = 0.0004 m=1.19, NMSE=0.27825 µ 0 = 0.49 m=2.715, NMSE=0.15338 µ 0 = 0.3, ν 0 = 0.05 m=2.75, NMSE=0.07523

Figure 8 .

 8 Figure 8. KAPA: Mean Squared Error with and without adaptation -Polynomial kernel (β = 2)

  The adapted element is the one with the least |α i,n |, β = 3 µ 0 = 0.3 m=2.315, NMSE=0.21308 µ 0 = 0.15, ν 0 = 0.00055 m=1.735, NMSE=0.21577 µ 0 = 0.3662 m=2.915, NMSE=0.15078 µ 0 = 0.3, ν 0 = 0.05 m=2.91, NMSE=0.12448

Figure 9 .

 9 Figure 9. KAPA: Mean Squared Error with and without adaptation -Polynomial kernel (β = 3)

Figure 10 .m= 9 .Figure 11 .

 10911 Figure 10. KNLMS: Mean Squared Error with and without adaptation -gaussian kernel

  µ 0 = 0.5 , m= 16.895 KRLS without adjustement, GRBF µ 0 = 0.5 , m= 16.895 KRLS with adjustement, GRBF µ 0 = 0.5, ν 0 = 0.2, m = 13.28 KNLMS with adjustement, GRBF µ 0 = 0.5, ν 0 = 0.085, m = 14.445

Figure 12 .Figure 13 .

 1213 Figure 12. KNLMS versus KRLS -gaussian kernel

Table I .

 I KAPA Experimental Setup and Performance, with η = 0.01, and ǫ = 0.07 -Polynomial kernel

		Third order polynomial kernel β = 2	Second order polynomial kernel β = 3
		Adaptation of the element with	Adaptation of the element with
		the highest norm of the gradient	the smallest value of |α n,i |
	Adaptation Parameters Settings	m	NMSE	Parameters Settings	m	NMSE
	Without	µ 0 =0.3	1.62 0.27649	µ 0 =0.3	2.315 0.21308
	Without	µ 0 =0.49	2.715 0.15338	µ 0 =0.3662	2.915 0.15078
	With	µ 0 =0.3, ν 0 =0.05	2.75 0.07523	µ 0 =0.3, ν 0 =0.05	2.91 0.12448
	With	µ 0 =0.15, ν 0 =0.0004 1.19 0.27825 µ 0 =0.15, ν 0 =0.00055 1.735 0.21577

Table II .

 II KNLMS Experimental Setup and Performance, with η = 0.09, and ǫ = 0.03

		Gaussian kernel		Exponential kernel
	Adaptation Parameters Settings	m	NMSE	Parameters Settings	m	NMSE
	Without	µ 0 =0.5	16.895 0.00628	µ 0 =0.3	9.495 0.019287
	Without	µ 0 =0.44	14.38	0.00845	µ 0 =0.31	10.255 0.012407
	With	µ 0 =0.5, ν 0 =0.1	14.36	0.0021	µ 0 =0.3, ν 0 =0.25	10.585 0.00324
	With	µ 0 =0.3, ν 0 =0.25	11.065 0.006262	µ 0 =0.2, ν 0 =0.012	6.5	0.01932

Table II gives a recapitulation for the obtained results and leads to the following observations: 23


Tian Wang
email: wangtian@buaa.edu.cn

Jie Chen
email: dr.jie.chen@ieee.org

Paul Honeine
email: paul.honeine@utt.fr

Hichem Snoussi
email: hichem.snoussi@utt.fr

Abnormal event detection via multi-kernel learning for distributed camera networks

Keywords: Distributed camera networks, abnormal detection, optical flow, one-class SVM

Abnormal event detection via multikernel learning for distributed camera networks

I. INTRODUCTION

Detecting abnormal events via video sequence analysis is crucial for public security management. In complex scenes, distributed camera networks with overlapping views are capable to obtain additional information to surveil the movement of the crowds compared to the single camera setting. In the PETS dataset [START_REF] Pets | Performance evaluation of tracking and surveillance (pets) 2009 benchmark data. multisensor sequences containing different crowd activities[END_REF], camera locations are illustrated in Fig. 1 directions, implying that people are attracted by some particular events, consequently, these scenes are considered as abnormal.

The scenes are captured by different cameras. Fig. 2(a)(b) are captured by camera 1 which is set at the side of the road, the movement of the people are well captured. Figs. 2(c)(d) are captured by camera 2 which faces the movement direction of the people, and there is occlusion of the individuals in this view. Figs. 2(e)(f) are captured by camera 3 which is also set at the side of the road, with larger distance. The purpose of the distributed camera surveillance is to detect abnormal events by benefitting the multi-view video sequences. Distribute camera networks are now widely used for surveillance application.

In the literature papers [START_REF] Wan | Distributed bayesian inference for consistent labeling of tracked objects in nonoverlapping camera networks[END_REF], [START_REF] Jiuqing | Multiple people tracking using camera networks with overlapping views[END_REF], the framework for multiple pedestrian tracking by using overlapping cameras was presented.

In [START_REF] Sankaranarayanan | Distributed sensing and processing for multi-camera networks[END_REF], several major challenges in distributed video processing, including robust and computationally efficient inference, opportunistic and parsimonious sensing were discussed. Large-scale video networks starts to play an important rule for video surveillance, object recognition, abnormal event detection and people tracking in crowded environments.

Modeling the movement feature of pixels is fundamental for detecting the abnormal event. In [START_REF] Cui | Abnormal detection using interaction energy potentials[END_REF], a method that tracked the local spatio-temporal interest points was proposed, and the abnormal activity was indicated by uncommon energy-velocity of the feature points. In [START_REF] Wang | Evaluation of local spatio-temporal features for action recognition[END_REF], a spatio-temporal descriptor was computed based on computing the histograms of optical flow in the neighborhood of detected points. In [START_REF] Wang | Action recognition by dense trajectories[END_REF], dense points were sampled from each frame and were tracked based on displacement information from an optical flow field. But for the crowd event analysis, it is difficult to obtain the pre-detected pixels of the blob due to the occlusion of the individuals.

In order to deal with the uncertainty of observations existing in video events, Bayesian modeling approaches such as hierarchical Dirichlet processes were used in [START_REF] Haines | Delta-dual hierarchical dirichlet processes: A pragmatic abnormal behaviour detector[END_REF], probabilistic latent semantic analysis was used in [9]. A method based on the variable-duration hidden Markov model was proposed in [10], where the durations of states were modeled except for the transitions between states, and the temporal understanding of the structure of complex events was tackled. Latent Dirichlet allocation (LDA) was also a typical standard topic model which has been used to model video clips as being derived from a bag of topics drawn from a fixed set of proportions [11]. In [12], the covariance matrix descriptor fusing the optical flow to encode moving information of a frame was presented. These work focused on the single view scene video analyzing, the feature abstraction method or the event models were heavily researched in the abnormal detection problem. Moreover, the abnormal event detection problems in distributed camera networks have a considerable room.

The rest of the paper is organized as follows. In Section 2, the optical flow-based feature is presented. In Section 3, the abnormal detection framework based on one-class SVM classification method is presented, and then, a multi-kernel strategy is proposed to deal with abnormal event detection problem for distributed camera networks. In Section 4, the experimental results are illustrated and discussed. Finally, Section 5 concludes the paper and gives a perspective of future work.

II. FEATURE SELECTION FOR ABNORMAL DETECTION

Horn-Schunck (HS) [START_REF] Horn | Determining optical flow[END_REF] is chosen to compute the optical flow which represents the movement information. The HS method formulates the optical flow as a global energy functional for the gray image sequence:

E = [(I x u + I y v + I t ) 2 + α( ∇u 2 + ∇v 2 )]dxdy (1)
where I x , I y and I t are the derivatives of the image intensity values along the horizontal direction x, vertical direction y and time t dimension, respectively. u, v are the horizontal and vertical optical flow. α is a regularization constant.

Based on the optical flow, the histogram of the optical flow orientation (HOFO) [START_REF] Wang | Detection of abnormal visual events via global optical flow orientation histogram[END_REF] is computed to fuse the movement as a

III. ONE-CLASS SVM WITH MULTIPLE KERNELS

The problem of non-linear one-class SVM [START_REF] Schölkopf | Estimating the support of a high-dimensional distribution[END_REF], [START_REF] Canu | Svm and kernel methods matlab toolbox[END_REF] can be cast as a quadratic programming problem:

min ω,ξ,ρ 1 2 w 2 + 1 νn n i=1 ξ i -ρ, subject to w, Φ(x i ) ≥ ρ -ξ i , ξ i ≥ 0. (2)
where x i ∈ X with i = 1 . . . n are n training samples in the original data space X . ξ i is the slack variable for penalizing the outliers. The hyperparameter ν ∈ (0, 1] is the weight for controlled slack variable, it tunes the number of acceptable outliers.

w, Φ(x i ) -b = 0 is the decision hyperplane. Φ is defined for building the non-linear classification problems, and it is a map from the set of the original input data X to a feature space H where the classification problem has a linear solution. The inner product in space H is defined by the kernel function κ(x i , x j ) = Φ(x i ) • Φ(x j ) . Introducing the Lagrangian multipliers α i , the decision function in the data space X is:

f (x) = sgn( n i=1 α i κ(x i , x) -ρ), (3)
where x is a vector in the input data space X , κ is the kernel function implicitly mapping the data into a high dimensional feature space where a linear classifier can be designed.

The Gaussian kernel is chosen to handle movement feature in this work. It is a semi-positive definite kernel that satisfies the Mercer condition [START_REF] Vapnik | The nature of statistical learning theory[END_REF], [START_REF] Vapnik | Statistical learning theory[END_REF], and defined by

κ(x i , x j ) = exp - x i -x j 2 2σ 2 , (x i , x j ) ∈ X × X , (4)
where x i , x j are the data in the original data space X , the variance σ indicates the scale factor at which the data should be clustered.

For constructing a more representative and discriminative feature descriptor for the distributed camera network, we take the scene captured by each view as a partial feature. The multi-kernel strategy considers the linear combination of candidates kernels:

κ(x i , x j ) = m s=1 µ s κ s (x i , x j ). (5)
where κ s , s = 1, • • • , m are m candidate kernels that satisfy the Mercer condition, and µ s are nonnegative factors. Consequently, their combination κ is also a semi-positive definite kernel. In this expression, the Gaussian kernel is adopted with:

κ s (x i , x j ) = exp - x i [s] -x j [s] 2 2σ 2 . ( 6
)
The kernels κ s , s = 1, • • • , m are the Gaussian kernels in this paper. Each sample vector

x consists of m parts [x 1 , x 2 , . . . , x m ].
For a given scene monitored by multi-camera, supposing that a set of training frames are obtained, based on the one-class SVM hypothesis, the abnormal behavior is the sample deviating from the training set. For example, a plaza is monitored by 3 cameras as show in Fig. 2. If s = 1, the frame captured by camera 1 is selected. We preset the weight µ s according to the characteristics of the image to tune the importance of each view. By considering this combination, the resulting kernel outperform each kernel κ s used individually. Based on the histogram of the optical flow orientation feature descriptor and the nonlinear one-class SVM, the abnormal event detection method is summarized in Algorithm 1, and explained as following.

Algorithm 1 Abnormal event detection algorithm.

Require:

Image set captured by the cameras.

1: Computing the optical flow of the training frame set

[I v i 1 , . . . , I v i m ],
where vi means the scene is monitored by camera i via the HS optical flow method:

[I v i 1 , . . . , I v i m ] -→ [O v i 1 , O v i 2 , . . . , O v i m ]
2: Computing the histogram of the optical flow orientation (HOFO) of the image in different views:

[O v i 1 , O v i 2 , . . . , O v i m ] -→ [H v i 1 , . . . , H v i m ] 3:
The feature sample of the image k under c distributed camera network is:

H k = [H v 1 k , H v 2 k , . . . , H vc k ] 4:
Training feature sample are learned by the nonlinear one-class SVM method to obtain the support vectors:

[H1, . . . , Hm] -→ support vector [S v i 1 , . . . , S v i m 1 ]
5: Each incoming frame Hp,...,q is classified by the decision function of one-class SVM.

6: The normal event or abnormal event is detected.

Step 1: The optical flow features are computed. The training frame set [I vi 1 , . . . , I vi m ] monitored by multi-camera network describing the normal behavior is available. I vi j means that the j-th image is captured by the i-th camera. The Horn-Schunck method is applied to acquire the optical flow feature. This step can be presented as:

[I vi 1 , . . . , I vi m ] -→ [O vi 1 , O vi 2 , . . . , O vi m ]
where Step 3: Fusing the 1-st to c-th HOFO feature to one high dimension feature vector. It is described in the equation below:

H k = [H v1 k , H v2 k , . . . , H vc k ]
where H k is the feature vector fusing the multi-view movement information. H v1 k is the HOFO feature of the k-th image in view i.

Step 4: Nonlinear one-class SVM is applied on the training frame HOFO descriptor in multi-view to obtain the support vector. It is described as the following:

[H 1 , . . . , H m ] -→ support vector [S 1 , . . . , S m1 ]
where Step 5: In the online detection phase, based on the support vectors obtained in the taring step, the one-class SVM classifies each incoming frame feature [H p , . . . , H q ]. Thus, the normal or abnormal event of the frame is classified, thus the abnormal event is detected.

IV. ABNORMAL EVENTS DETECTION RESULTS

We then conduct experiments to evaluate the performance of the one-class SVM classification method for abnormal frame event detection with a distributed camera network. The PETS [START_REF] Pets | Performance evaluation of tracking and surveillance (pets) 2009 benchmark data. multisensor sequences containing different crowd activities[END_REF] dataset is used for the evaluation purpose. In the experiments of PETS [START_REF] Pets | Performance evaluation of tracking and surveillance (pets) 2009 benchmark data. multisensor sequences containing different crowd activities[END_REF] dataset, each event is represented by 3 separately scenes, thus, the event is described by 3 separated HOFO features which is marked as "3 HOFO". Otherwise, we mark the feature as "1 HOFO". If the multi-kernel strategy is used, we mark it as "3 kernels", otherwise, we mark the kernel strategy as "1 kernel". The detection accuracy of the detection results are shown for the experiments.

The normal and abnormal event of sequence Time 14-17 in the PETS dataset are shown in Fig. 4. The training samples and the normal samples for testing are the frames in 3 views chosen from the sequence (Time 14-55) where the people are walking in different directions. 400 training frames (Frame0000 to Frame0399), and 90 normal testing frames (Frame0400 to Frame0489) are selected from Time14-55. The abnormal testing samples are selected from the sequence (Time 14-17) where the people are walking or running in the same direction. 89 abnormal testing frames (Frame0000 to Frame0089) are selected from Time14-17.

The AUC (area under the curve) value of the abnormal detection results in different views and different multi-kernel strategies are shown in TABLE ??. "Single View" means "1 HOFO 1 kernel" strategy, and "Multi-view" means "3 HOFO 3 kernels".

"View 1" means the scene is monitored by camera 1, the abnormal detection results are shown in Figs. View 2 has the lowest area under the ROC (receiver operating characteristic curve), since as previously mentioned, the camera 2 faces the movement direction of the crowd, the occlusion influences the computation of the optical flow. Thus the HOFO feature based on the optical flow cannot represent the accurate movement information. The results show that the abnormal detection algorithm of HOFO feature can obtain satisfactory detection results. Moreover, the multi-kernel strategy can generally improve the performance. The parameters of the multi-kernel influence the performance. Thus, the optimal parameters should be sought.

The experiment detecting the running activity as the abnormal event is shown in Fig. 6 and Fig. In the future work, the optimal coefficients of the multi-kernel strategy should be obtained automatically based on the scene, while these parameters were pre-set in current work. Additional datasets will be considered to show the advantages of distributed camera networks, such as single person action recognition or action tracking in multi-view scenes, etc.

detect:walk

  . Several normal and abnormal scenes from different angles of view are shown in Fig. 2. In Figs. 2(a)(c)(e), all people are walking in different directions, which are considered as normal. In Figs. 2(b)(d)(f), people are walking or running towards the same

Fig. 1 :

 1 Fig. 1: The plan of the multi-camera localizations in the PETS dataset. Three cameras are set in the campus.

Fig. 2 :

 2 Fig. 2: Examples of the normal and abnormal scenes captured by distributed cameras of PETS dataset. (a)(c)(e): The people are walking in different directions, the normal scenes of Time 14-55 sequence; (b)(d)(f): The people are walking in the same direction, the abnormal scenes of Time 14-17 sequence; (a)(b): scenes captured by camera 1; (c)(d): scenes captured by camera 2; (e)(f): scenes captured by camera 3.

  high dimension vector. A 2 × 2 rectangular cell HOFO descriptor of the image is shown in Fig. 3. The orientation is computed by horizontal and vertical optical flow, and then it is voted into n bins. Each cell contains h c pixels in height dimension, and w c in width dimension. A block contains h b × w b cells, it is set as 2 × 2 in this paper.

Fig. 3 :

 3 Fig. 3: Histogram of the optical flow orientation (HOFO) feature descriptor of the image.

  [H 1 , . . . , H m ] are the histogram of optical flow orientation descriptors of the training frames under multi-view environment. [S 1 , . . . , S m1 ] are support vectors that are the minority of the training vectors contribute to the decision function.

Fig. 4 :

 4 Fig. 4: Detection results of the normal and abnormal scenes of Time 14-17 captured by distributed camera network. (a)(c)(e): The detection result of one normal frame. (b)(d)(f): The detection result of one abnormal frame.

7

 7

Fig. 5 :

 5 Fig. 5: ROC curve of abnormal detection result of sequence Time 14-17 under different views or different kernel strategies.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Detection results of the normal and abnormal scenes of Time 14-16 captured by thedistributed camera network, the individuals are moving from right to left. (a)(c)(e): The detection result of one normal frame, the individuals are walking in one direction. (b)(d)(f): The detection result of one abnormal frame, the individuals are running in one direction.

TABLE I :

 I The abnormal detection results of sequence Time 14-17. The comparison of the abnormal frame event detection results in single view scene and in distributed camera scenes via multi-kernel strategies. in PETS dataset where the individuals are walking in one direction. In this experiment, 61 frames (Frame0000 to Frame0060) where people are walking from left to right in sequence Time14-17, and 50 frames (Frame0000 to Frame0049) where the individuals are walking from right to left in sequence Time14-31. Corresponding, 104 normal samples and 118 abnormal frames in sequence Time 14-16 are detected. The abnormal detection performance is improved by multi-kernel strategy also. For example, the AUC values of the the sequence where the individuals are moving from left to right are shown

	Strategy	Area under ROC
		View 1 View 2	View 3
	Single View	0.9036	0.7807	0.9270
	Multi-view(1)		0.9441
	Multi-view(2)		0.9643
	Time14-31)

in TABLE ??.

TABLE II :

 II The abnormal detection results of sequence Time 14-16. The comparison of the abnormal frame event detection results in single view scene and in distributed camera scenes via multi-kernel strategies. abnormal frame event detection with distributed camera networks is proposed in this paper. The histogram of optical flow is computed as the descriptor to represent the movement of a frame. A multi-kernel strategy is presented to benefit the characteristics of distributed camera network, and the performance of the abnormal detection results via nonlinear one-class SVM is improved. The benchmark dataset PETS has been tested to demonstrate the effectiveness of the proposed algorithm.

	Strategy	Area under ROC
		View 1 View 2	View 3
	Single View	0.8637	0.8071	0.9205
	Multi-view(1)		0.9312
	Multi-view(2)		0.9403
	V. CONCLUSIONS
	A method for

ACKNOWLEDGMENT

This work is partially supported by the SURECAP CPER project (fonction de surveillance dans les réseaux de capteurs sans fil via contrat de plan Etat-Région) and the Platform CAPSEC (capteurs pour la sécurité) funded by Région Champagne-Ardenne and FEDER (fonds européen de développement régional), the Fundamental Research Funds for the Central Universities and the National Natural Science Foundation of China (Grant No. U1435220).

Conference on Computer Vision (ICCV), 2011, pp. 2198-2205.

[9] J. Varadarajan and J.-M. Odobez, "Topic models for scene analysis and abnormality detection," in Proceedings of the 12th International Conference on Computer Vision Workshops (ICCV Workshops), 2009, pp. 1338-1345.

[10] K. Tang, L. Fei-Fei, and D. Koller, "Learning latent temporal structure for complex event detection," in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (CVPR). IEEE, 2012, pp. 1250-1257.

[11] O. P. Popoola and K. Wang, "Video-based abnormal human behavior recognition-a review," IEEE Transactions on Systems, Man, and Cybernetics, Part C: Applications and Reviews, vol. 42, no. 6, pp. 865-878, 2012.

[12] T. Wang, J. Chen, Y. Zhou, and H. Snoussi, "Online least squares one-class support vector machines based abnormal visual event detection," Sensors, no. 12, pp. 17 130-17 155, 2013.