Latifa Oukhellou

Etienne Côme

Patrice Aknin

Thierry Denoeux

Semi-supervised feature extraction using independent factor analysis

Keywords: Independent factor analysis, semi-supervised learning, mixture models, maximum likelihood, dimensionality reduction

come

Introduction

A considerable amount of research has been devoted to dimensionality reduction of multivariate data sets. The underlying motivation is that many data sets live on a subspace whose intrinsic dimensionality is lower than that of the original data space. Research in this area has employed either feature selection methods which directly select a subset of meaningful variables from the original ones, or feature extraction methods that aim to generate new features from the first input data. Whatever the chosen approach, the goal is to describe a large data set by a reduced number of variables that better capture the essential structure of the problem. These methods include linear approaches such as probabilistic principal component analysis (PPCA) [START_REF] Jolliffe | Principal Component Analysis[END_REF][START_REF] Jackson | A User's Guide to Principal Components[END_REF] Projection Pursuit (PP) [START_REF] Friedman | Exploratory projection pursuit[END_REF][START_REF] Jones | What is projection pursuit ?[END_REF], Metric Muldimensional Scaling (MDS) [START_REF] Cox | Multidimensional Scaling[END_REF], linear discriminant analysis (LDA) [START_REF] Fukunaga | Introduction to statistical pattern recognition[END_REF] on the one hand, and non linear methods such as random projection (RP) [START_REF] Johson | Extensions of Lipschits Mapping into a Hilbert Space[END_REF][START_REF] Watanabe | Dimensionality Reduction by Random Projection[END_REF], Kohonen's self-organizing maps [START_REF] Kohonen | The self-organizing map[END_REF][START_REF] Kaski | Data exploration using self-organizing maps[END_REF],... on the other hand. A detailed survey of many of these methods can be found in [START_REF] Saul | Spectral methods for dimensionality reduction[END_REF][START_REF] Burges | Geometric methods for feature extraction and dimensional reduction[END_REF].

The dimensionality reduction problem can also be formulated using a generative latent variable model which aims to describe observed variables (original features), in terms of smaller set of unobservable (or latent) variables. Depending on the assumption made on the latent and observed variable Email address: oukhellou@inrets.fr Tel: +33 1 45 92 56 58 (L. Oukhellou) distributions, different kind of models can be distinguished such as Principal Component Analysis (PCA), Factor analysis (FA) [START_REF] Bartholomew | Latent variable models and factor analysis. Kendall's library of satistics[END_REF][START_REF] Everitt | An Introduction to Latent Variable Models[END_REF], and Independent Component Analysis (ICA) [START_REF] Hyvärinen | Independent Component Analysis[END_REF][START_REF] Bell | An information maximization approach to blind separation and blind deconvolution[END_REF][START_REF] Amari | A New Learning Algorithm for Blind Signal Separation[END_REF]. ICA has been applied to many different problems, including blind source separation, exploratory data analysis and feature extraction [START_REF] Cichocki | Adaptive Blind Signal and Image Processing[END_REF][START_REF] Bakir | Non-negative maximum likelihood ICA for blind source separation of images and signals with application to hyperspectral image subpixel demixing[END_REF]. In the feature extraction context, several authors used ICA to extract meaningful features for both regression and classification problems [START_REF] Kwak | Dimensionality reduction based on ICA for regression problems[END_REF][START_REF] Kwak | Feature extraction based on ICA for binary classification problems[END_REF]. This paper deals with a particular model of this family, recently proposed by [START_REF] Attias | Independent factor analysis[END_REF][START_REF] Moulines | Maximum likelihood for blind separation and deconvolution of noisy signals using mixture models[END_REF], and known as Independent Factor Analysis (IFA).

The generative model involved in IFA assumes that observed variables are generated by a linear mixture of independent and non Gaussian latent variables as in the ICA model. Furthermore, it considers that each individual latent variable has its own distribution, modeled by a mixture of Gaussians (MOG). The IFA model is often considered within an unsupervised learning framework. The model parameters and thus the latent variables are learned from the observed data only. Recent works have derived an approach for modeling class conditional densities based on IFA model [START_REF] Montanari | Independent factor discriminant analysis[END_REF]. In this paper, we propose an extension of the basic IFA model makes it possible to incorporate additional information on cluster membership of some training samples to estimate the IFA model. In this way, the learning of this model, and thus the dimensionality reduction can be handled in a semi-supervised learning framework.

The article is organized as follows. In Section 2, we review the independent factor analysis model and present how it can be estimated by maximum likelihood in a noiseless setting. Section 3 focuses on the problem of semi-supervised learning of the IFA model where additional information on cluster membership of some samples will be incorporated. A generalized maximum likelihood criterion will be defined and the algorithm for its optimization also detailed. Experimental results showing the benefits of the proposed approach to achieve dimensionality reduction will then be given for real data sets. Conclusions are presented in Section 5.

Noiseless Independent Factor Analysis

Background on Independent Factor Analysis

IFA was introduced in [START_REF] Moulines | Maximum likelihood for blind separation and deconvolution of noisy signals using mixture models[END_REF][START_REF] Attias | Independent factor analysis[END_REF]. It originates from both ordinary factor analysis (FA) in applied statistics [START_REF] Spearman | General intelligence, objectively determined and measured[END_REF][START_REF] Thusrstone | Multiple Factor Analysis[END_REF] and Independent Component Analysis (ICA) in signal processing [START_REF] Hyvärinen | Independent Component Analysis[END_REF][START_REF] Bell | An information maximization approach to blind separation and blind deconvolution[END_REF]. IFA aims to recover independent latent variables from their observed linear mixtures. The latent variables are assumed to be mutually independent and non Gaussians. In the noiseless form that is used throughout this paper, the IFA model can be expressed as:

x = A z, (1)
where A is a square matrix of size S × S , x the random vector whose elements (x 1 , . . . , x S ) are the mixtures and z the random vector whose elements (z 1 , . . . , z S ) are the latent variables. Thanks to the noiseless setting, a deterministic relationship between the distributions of observed and latent variables can be expressed as:

f X (x) = 1 | det(A)| f Z (A -1 x), (2)
Unlike the ICA model in which the probability density functions of the latent variables are fixed using prior knowledge or according to some indicator that allows switching between sub and super Gaussian densities [START_REF] Hyvärinen | Independent Component Analysis[END_REF], each latent variable density in the IFA is modeled as a mixture of normally distributed components (Mixture of Gaussians MOG) so that a wide class of densities can be approximated [START_REF] Moulines | Maximum likelihood for blind separation and deconvolution of noisy signals using mixture models[END_REF][START_REF] Attias | Independent factor analysis[END_REF][START_REF] Côme | Apprentissage de modèles génératifs pour le diagnostic de systèmes complexes avec labellisation douce et contraintes spatiales[END_REF]:

f Z s (z s ) = K s k=1 π s k ϕ(z s ; µ s k , ν s k ), (3)
where ϕ(.; µ, ν) denotes a univariate normal density function with mean µ and variance ν. Equation (3) means that each latent variable is described as mixture of K s Gaussians with mean µ s k , variance ν s k and mixing proportions π s k . Considering the graphical model of IFA shown in Figure 1, it can be seen that the IFA model provides two levels of interpretation corresponding to discrete and continuous latent variables. For each one of the S latent variables, the discrete latent variable encodes the cluster from which each sample is drawn.

The whole IFA model parameters can thus be summarized in a vector ψ = (W, π 1 , . . . , π S , µ 1 , . . . , µ S , ν 1 , . . . , ν S ), with W the unmixing matrix (W = A -1 ), π j the vector of cluster proportions of source j which sum to 1, µ j and ν j the vectors of size K j containing the means and the variances of each cluster.

Continuous latent variables

Discrete latent variables

Figure 1: Graphical model for Independent Factor Analysis.

Parameter estimation in IFA 2.2.1. IFA and Maximum Likelihood

The learning problem associated with the IFA model consists in estimating both the unmixing matrix W and the MOG parameters from the observed variables alone. Considering an i.i.d random sample X = (x 1 , ..., x N ) of size N and using Equation (2) under the latent variable independence hypothesis, the log-likelihood has the form:

L(ψ; X) = N i=1 S s=1 log f Z s ((Wx i ) s ) + N log(| det(W)|), (4)
By substituting the density distribution by its expression given in (3), the log-likelihood can be rewritten as:

L(ψ; X) = N log(| det(W)|)+ N i=1 S s=1 log         K s k=1 π s k ϕ (W -1 x i ) s , µ s k , ν s k         .
(5) The whole IFA parameters ψ can therefore be estimated by maximizing the likelihood function.

Generalized Expectation-Maximization (GEM) Algorithm

When the latent variable densities are known, the unmixing matrix estimation is based on gradient methods that maximize the likelihood. The gradient of the log-likelihood defined in ( 5) can be derived as:

∂L(W; X) ∂W ∝ (W -1 ) t - 1 N N i=1 g(Wx i )x t i , (6)
where

g(z) = -∂ log( f Z 1 (z 1 )) ∂z 1 , . . . , -∂ log( f Z S (z S )) ∂z S t . (7)
The update rule of the unmixing matrix is thus given by:

W (q+1) = W (q) + τ        (W (q) ) -1 t - 1 N N i=1 g(W (q) x i )x t i        , ( 8
)
where τ is the gradient step that can be adjusted by means of linear search methods ( [START_REF] Nocedal | Numerical Optimization[END_REF]).

The convergence of this algorithm can be improved by using the natural gradient ([16, p. 67, p. 208]).

Maximum likelihood of the whole model parameters can be achieved by an alternating optimization strategy. The gradient algorithm is indeed well suited to optimize the log-likelihood function with respect to the unmixing matrix W when the parameters of the source marginal densities are frozen. Conversely, with W kept fixed, an EM algorithm can be used to optimize the likelihood function with respect to the parameters of each source. These remarks have led use a Generalized EM algorithm (GEM) [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF][START_REF] Mclachlan | The EM algorithm and Extension[END_REF] that simultaneously maximizes the likelihood function with respect to all model parameters.

3. Semi-supervised learning in Independent Factor Analysis

Derivation of a Generalized Likelihood Criterion

The IFA model is often considered within an unsupervised learning framework. This section considers the learning of this model in a partially-supervised learning context where partial knowledge of the cluster membership of some samples is available. For that purpose, the model is built from a combination of M labeled and N -M unlabeled samples. Consequently, the criterion can be decomposed into two parts corresponding, respectively, to the supervised and unsupervised learning examples and the log-likelihood criterion (5) can be written as:

L(W; X) = N log(| det(W|)+ M i=1 S s=1 K s k=1 l s ik log π s k ϕ (Wx i ) s , µ s k , ν s k + N i=M+1 S s=1 log         K s k=1 π s k ϕ (Wx i ) s , µ s k , ν s k         . ( 9
)
Note that l s ik ∈ {0, 1} K s , l s ik = 1 if sample i comes from component c k of sources s and l s ik = 0 otherwise.

Practical Considerations

A Generalized EM algorithm (GEM), (also noted here as Algorithm 1) can be designed to simultaneously maximize the likelihood function with respect to all the model parameters [START_REF] Côme | Partiallysupervised learning in Independent Factor Analysis In 17 th[END_REF][START_REF] Côme | Noiseless Independent Factor Analysis with mixing constraints in a semi-supervised framework. Application to railway device fault diagnosis[END_REF]. This algorithm is similar to the EM algorithm used to estimate IFA parameter in an unsupervised setting [START_REF] Attias | Independent factor analysis[END_REF], except for the E step, where the posterior probabilities t s ik are only computed for the unlabeled samples. The score function g of each latent variable density are given by:

g s (z is ) =          -∂ log Ks k=1 l s ik π s k ϕ(z is ;µ s k ,ν s k ) ∂z is , if i ≤ M -∂ log Ks k=1 π s k ϕ(z is ;µ s k ,ν s k ) ∂z is , if i > M =          K s k=1 l s ik (z is -µ s k ) ν s k , if i ≤ M K s k=1 t s ik (z is -µ s k ) ν s k , if i > M ( 10
)
Note that t s ik is the posterior probability that the sample i belongs to component k of the latent variable s, given z is = (Wx i ) s and the labels:

t s ik = π s k ϕ(z is ; µ s k , ν s k ) K s k =1 π s k ϕ(z is ; µ s k , ν s k ) . ( 11
)

Simulations and discussion

In this section, we investigated the interest of our approach with three real datasets of which the characteristics are given on Table 1. The first dataset is the Crabs dataset1 that concerns the recognition of crabs species and sexes in a population of crabs using different morphological measurements. The second dataset is the well known Fisher's Iris data available on-line 2 . The third dataset is the YaleB face database [START_REF] Georghiades | From Few to Many: Illumination Cone Models for Face Recognition under Variable Lighting and Pose[END_REF]. We choose the first 5 subjects from the dataset and get totally 320 face samples that were captured under different illumination conditions (5 subjects × 64 illumination conditions).

To better understand our approach as compared to the unsupervised IFA model, different experiments were carried out to show the influence of learning the IFA model when information regarding the component membership of some training samples is introduced. We show that such information can be exploited to efficiently extract a reduced set of variables from the original ones. We show also the potential benefit of incorporating labels in terms of simplification of the optimization problem. This consideration has high practical interest as the problem of local maxima is very important for independent factor analysis.

Three different learning strategies were compared, namely, PCA, unsupervised IFA, and semi-supervised IFA using some labelled samples over all latent variables. The IFA model provides two levels of interpretation corresponding to discrete and continuous latent variables. While results for continuous variables are useful to visualize the projection of the data onto the two-dimensional principal subspace, those of discrete variables allow quantification of classification rates according to each latent variable. In this case, a Maximum a posteriori (MAP) criterion is used to assign each observation to one of the mixture components modeling each latent variable density function. When PCA is applied to reduce dimensionality, the nearest neighborhood (1-NN) classifier is employed for classification.

A few data samples from the whole data are randomly chosen as labeled samples. Twenty random starting points were used for the GEM algorithm and only the best solution according to the likelihood was kept. The performances were quantified using the correct detection rates according to each latent variable calculated on a test set constituted by the samples not labelled during the training phase. The process is repeated for 30 runs and the averaged results are recorded.

Algorithm 1: Pseudo-code for Semi-supervised IFA with GEM algorithm Input: Centered observation matrix X, cluster membership for the M labeled data l s ik # Random initialization of IFA parameter vector ψ (0) , q = 0 while convergence test do # latent variable update Z = X.W (q) t # Update of the latent variable parameters / EM forall s ∈ {1, . . . , S } and k ∈ {1, . . . ,

K s } do # E-Step t s(q) ik = π s(q) k ϕ(z is ; µ s(q) k , ν s(q) k ) K s k =1 π s(q) k ϕ(z is ; µ s(q) k , ν s(q) k ) , ∀i ∈ {1, . . . , N -M} t s(q) ik = l s ik , ∀i ∈ {1, . . . , M} forall s ∈ {1, . . . , S } and k ∈ {1, . . . , K s } do # M-step , Update of the parameter vector of each latent variable π s(q+1) k = 1 N N i=1 t s(q) ik µ s(q+1) k = 1 N i=1 t s(q) ik N i=1 t s(q) ik z is ν s(q+1) k = 1 N i=1 t s(q) ik N i=1 t s(q) ik (z is -µ s(q+1) k ) 2
# Update of the score matrix G (10)

G = g (q+1) (Z) # Gradient (8) ∆W = ( (W (q) ) -1 t -1 N N i=1 g(W (q) x i )x t i # Linear search τ (gradient step) τ * = Linearsearch(W (q) ,∆W) # Unmixing matrix update W (q+1) = W (q) + τ * .∆W # Latent variable normalization to remove scale indetermination forall s ∈ {1, . . . , S } do σ 2 s = K s k=1 π s(q+1) k (ν s(q+1) k +µ s(q+1) k 2 )-K s k=1 π s(q+1) k µ s(q+1) k 2 forall k ∈ {1, . . . , K s } do µ s(q+1) k = µ s(q+1) k /σ s ν s(q+1) k = ν s(q+1) k /σ 2 s W (q+1) s. = W (q+1) s.
/σ s q ← q + 1 Output: Estimated parameters : ψ, estimated latent variables : Z

Crabs dataset

The crabs dataset consists of 5 morphological measurements recorded for 200 crabs that can be categorized into four groups on the basis of their sex and species: "Blue male", "Blue female", "Orange Male", and "Orange female". The IFA model used to deal with this dataset has 5 latent variables, among which two were modelled by a mixture of 2 normally distribution components and were labelled by using the sex and the species information, the remaining latent variables are modeled by simple Gaussians.

Figure 2 shows the scatter plot of the crabs data when they are projected onto the estimated two-dimensional principal subspace obtained by the semi-supervised IFA model using 20% of labelled training data. As a comparison, the projection of the data onto the first two principal components (PCA) is also given. The semi-supervised IFA leads to a projection giving much better class separation. In fact, the first latent variable clearly captures the variability of the crabs sex while the second corresponds to their species.

First latent variable

Semi-supervised IFA

Second latent variable

Figure 2: Two-dimensional visualization of the Crabs dataset projected onto the first two principal components obtained by PCA and semi-supervised IFA with 20% of labelled samples. In the graph, the stars denote the "blue male" group of crabs, the circles denote the "blue female" group, the boxes indicate the "orange male" crabs and the triangles the "orange female".

In order to accurately quantify the affect of the proportion of labelled samples, the correct detection rates (for the 2 latent variables) has been evaluated as a function of the number of labelled samples. Figure 3 shows the classification performance when the proportion of labelled samples increases from 0 (unsupervised learning case) to 80%. With only 20% of labelled samples, the correct detection rates reach 91.4% for the sex latent variable and 100% for the species latent variable. Note that unsupervised IFA classify correctly classifies 42% of observations for the sex and 56% for the species. PCA-1NN gives 94.8% and 51.8% of correct classification rates for the sex and species. The results are summarized on Table 2.

Iris dataset

The Iris dataset consists of 50 samples from each of three species of Iris flowers (Iris setosa, Iris virginica and Iris versicolor). 4 features were measured from each sample : the length and the width of sepal and petal. The IFA model used for this dataset has only one latent variable with a mixture density of three components, (one for each species), the remaining variables being as usual simple Gaussians.

Figure 5 shows the 2D data visualization obtained by the PCA and the semi-supervised IFA with 10% of labelled training samples. It can be seen that semi-supervised IFA requires one latent variable with three components to highlight the latent structure of the Iris dataset. The less separable classes are the virginica and the versicolor species while the setosa is the best predictible class. In terms of classification, 86% of samples are well classified in the principal subspace given by the PCA-1NN while only 33% of correct detection rate achieved by unsupervised IFA. This is mainly due to latent variable permutation. By permuting the mixture components a correct detection rate of 97.1% is reached by unsupervised IFA. This finding suggests that the species structure is very clear in this dataset.

-4 -3 -2 -1 0 1 2 3 4 -3 -2 -1 0 1 2 3 First latent variable PCA -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 3 4
First latent variable Semi-supervised IFA

Second latent variable

Figure 5: Two-dimensional visualization of the Iris dataset projected onto the first two principal components obtained by PCA and semi-supervised IFA with 10% of labelled samples. In the graph, the stars denote the "setosa" group, the circles denote the "versicolor" group and the triangles the virginica" species.

Figure 6 shows both the correct detection rate obtained according to the first latent variable and the CPU time computation as functions of the proportion of labelled samples. It can be seen that the performances are drastically improved when 10% of samples are labelled and slows down with more labelled data. This means that too many labelled data are not useful to improve the results for the Iris dataset because as already noted the species structure is very clear in this dataset. The results are also summarized in Table 2.

YaleB face dataset

5 subjects have been chosen from the YaleB dataset. For each one of them, 64 face samples have been captured under different illumination conditions. The complexity of dimensionality reduction task is due to the lightning conditions that lead to variability between images of the same subject greater than that of different subjects. For this dataset a pre-processing step was necessary to use semi-supervised IFA. In fact considering the pixels gray level as the input variable of the model will lead to a model with 32256 variables one for each pixel. We therefore first process the data by PCA and kept only the 40 leading principal components. The IFA models was then fixed as follows : 5 latent variables with mixture densities and 35 Gaussian latent variables. Each mixture has 2 components; one component encoding the membership of a picture to a specific subject and the remaining components encoding the membership of the picture to the remaining subjects. Therefore each latent variable with a mixture density can be used to classify picture according to one subject against all the others. Such modeling is partly in contradiction with the hypothesis of Independent Factor Analysis but we will see that good results are still obtained.

When the data are projected onto the first principal subspace given by the PCA, an important overlapping between the 5 groups is noticed. This is illustrated in Figure 7. However, it can be seen on the same figure that semi-supervised IFA with 40% of labelled samples captures much better the intrinsic structure of the dataset. By projecting the dataset on the latent variables corresponding to subjects 1 and 2 one can see that these variables are very good to distinguish these subjects from the others. The complete results averaged on the 5 subjects are summarized on Table 2. Figure 8 shows the boxplot of correct detection rate obtained for each subject when the proportion of labelled samples increases from 0 (unsupervised learning case) to 80%. With 40% of labelled samples, the classification performance reach 87% while it is equal to 50% in the unsupervised IFA.

Conclusion

In this paper, we have proposed a new approach of dimensionality reduction based on partially supervised Independent Factor Analysis. We introduced a criterion where the likelihood of both labelled and unlabelled samples is maximized by a GEM algorithm. Experimental results show that efficient dimensionality reduction can be achieved for some problems where unsupervised methods fail to capture the underlying structure of the data. The amount of labelled data required to reach a satisfactory accuracy depends on the discrimination problem complexity but for all the problems investigated here we find that with only 10% of labelled data semi-supervised IFA gives results much better than PCA and ordinary IFA.

A. Appendix : Gradient of the Log-likelihood with respect to the unmixing matrix

The log-likelihood in the noiseless IFA is given by:

L(ψ; X) = N i=1 S s=1 log f Z s ((Wx i ) s ) + N log(| det(W)|),
In order to compute the gradient of L(ψ; X) with respect to W, we have to compute the derivative of the logarithm of the absolute value of a matrix determinant with respect to one of its elements, which is given by (see [START_REF] Mackay | Maximum likelihood and covariant algorithms for independent component analysis[END_REF], [3, p. 8]):

∂ log(| det(X)|) ∂X lk = (X -1 ) kl , By using this relationship and assuming that latent variable densities f Z 1 , . . . , f Z S are known, the derivative of the loglikelihood with respect to one element W lk of the unmixing matrix can be written as:

∂L(ψ; X) ∂W lk = N(W -1 ) kl + N i=1 ∂ log f Z l ((Wx i ) l ) ∂W lk = N(W -1 ) kl - N i=1
x ik g l ((Wx i ) l ) ,

where g l (z) is the opposite of the derivative of the logarithm of the latent variable density z l : g l (z) = -∂ log( f Z l (z)) ∂z By using matricial notations, we can define the g function:

g : R S → R S g(z) = -∂ log( f Z 1 (z 1 )) ∂z 1 , . . . , -∂ log( f Z S (z S )) ∂z S t .

Which allows us to obtain the matrix of the derivative of the log-likelihood with respect to each element of W :

∂L(ψ; X) ∂W = N(W -1 ) t - N i=1 g (Wx i ) x t i ∝ (W -1 ) t - 1 N N i=1 g (Wx i ) x t i .

8 Figure 3 :

 83 Figure 3: Influence of the proportion of labelled samples on the estimation of the semi-supervised IFA model: Boxplot of the correct detection rates for the first (crab sex) and the second latent (crab species) variables function of the percentage of labelled samples.

Figure 4 Figure 4 :

 44 Figure 4 displays the CPU time required for the GEM algorithm convergence function of the amount of labelled samples.It can be seen that the time computation (or the number of iterations) exponentially decreases when the amount of labelled samples increases. This graph highlights the potential benefit of incorporating labels in terms of simplification of the optimization problem.

Figure 6 :

 6 Figure 6: Influence of the proportion of labelled samples on the estimation of the semi-supervised IFA: (a) Boxplot of the correct detection rates evaluated with the MAP criterion on the first latent variable. (b) Time computation over 20 random initializations for the GEM algorithm as a function of the proportion of labelled data

Figure 7 :

 7 Figure 7: Two-dimensional visualization of the YaleB face dataset projected onto the first two principal components obtained by PCA and semi-supervised IFA with 40% of labelled samples. In the graph, different symbols are used to denote the 5 subjects.

8 Figure 8 :

 88 Figure 8: Influence of the proportion of labelled samples on the estimation of the semi-supervised IFA: (a) Boxplot of the correct detection rates evaluated with the MAP criterion on the five latent variables. (b) Mean time computation over 20 random initializations for the GEM algorithm as a function of the proportion of labelled data

Table 1 :

 1 Characteristics of real datasets.

	name # dimensions # samples # classes
	Crabs	5	200	4
	Iris	4	150	3
	YaleB	40	320	5

Table 2 :

 2 Correct detection rates (%) averaged over 30 independent datasets, for three dimensionality reduction methods PCA-1NN, IFA and semi-supervised IFA with different proportions of labelled samples (%).

		PCA-1NN IFA			Semi-supervised IFA
		0	0	10	20	30	40	50	60	70
	Crabs	55.8	49.5 94.7 95.7 95.4 95.8	96	95.8 96.2
	Iris	86.1	33.1 97.4 97.1 97.5 97.3 97.1 97.1 97.1
	YaleB	25.8	49.7	71	75.7 80.1 84.1 86.9 87.5	89

http://rweb.stat.umn.edu/R/library/MASS/html/crabs.html

http://mlearn.ics.uci.edu/MLRepository.html


M Acheritogaray
email: marion.acheritogaray@math.univ-toulouse.fr

P Degond

A Frouvelle
email: amic.frouvelle@math.univ-toulouse.fr

J-G Liu
email: jliu@phy.duke.edu

Kinetic formulation and global existence for the Hall-Magneto-hydrodynamics system

Keywords: Hall-MHD, kinetic formulation, entropy dissipation, generalized Ohm's law, incompressible viscous flow, resistivity, global weak solutions, KMC waves AMS Subject classification: 35L60, 35K55, 35Q80

This paper deals with the derivation and analysis of the the Hall Magneto-Hydrodynamic equations. We first provide a derivation of this system from a twofluids Euler-Maxwell system for electrons and ions, through a set of scaling limits. We also propose a kinetic formulation for the Hall-MHD equations which contains as fluid closure different variants of the Hall-MHD model. Then, we prove the existence of global weak solutions for the incompressible viscous resistive Hall-MHD model. We use the particular structure of the Hall term which has zero contribution to the energy identity. Finally, we discuss particular solutions in the form of axisymmetric purely swirling magnetic fields and propose some regularization of the Hall equation.

Introduction

This paper deals with the derivation and analysis of the the Hall Magneto-Hydrodynamic (Hall-MHD) equations. The Hall-MHD model currently receives an increasing attention from plasma physicists. It is believed to be the key for understanding the problem of magnetic reconnection. Indeed, space plasma observations provide strong evidence for the existence of frequent and fast changes in the topology of magnetic field lines, associated to violent events such as solar flares [START_REF] Forbes | Magnetic reconnection in solar flares[END_REF]. However, magnetic reconnection cannot be described in the framework of ideal MHD, due to the frozen-field effect. Indeed, in ideal MHD, due to the Faraday equation and ideal Ohm's law, the magnetic field is essentially passively transported by the fluid velocity. Therefore, the topology of the magnetic field is preserved, even in the magnetic field lines are deformed by the flow. In order to break this passive magnetic field transport by the fluid flow, one is led to re-introduce the Hall terms which was neglected in ideal MHD. In spite of its increasing importance for physical applications, the Hall-MHD model has received very little attention from the theoretical viewpoint (see e.g. [START_REF] Mahajan | Exact solution of the incompressible Hall magnetohydrodynamics[END_REF][START_REF] Simakov | Quantitative, comprehensive, analytical model for magnetic reconnection in Hall magnetohydrodynamics[END_REF]) and the purpose of this paper is mainly to propose a framework for the derivation and analysis of the Hall-MHD problem.

We first provide a derivation of this system from a two-fluids isothermal Euler-Maxwell system for electrons and ions, through a set of scaling limits. The two-fluids model of plasma is known as Braginskii model [START_REF] Braginskii | Transport processes in a plasma[END_REF] and was justified in [START_REF] Duvaut | inéquations en thermoélasticité et magnétohydrodynamique[END_REF] (see also [START_REF] Degond | Asymptotic continuum models for plasmas and disparate mass gaseous binary mixtures[END_REF]) on the basis of a concurrent hydrodynamic and zero electron to ion mass ratio limit. As usual in MHD models, zero electron to ion mass ratio limit, zero Debye length limit and zero displacement current limits have to be taken. Then, the main point is to examine the orders of magnitude of the various terms arising in the generalized Ohm law (which is the electron momentum equation with zero inertia) and in the current equation. Basically, Hall MHD is obtained when the electron and ion velocities have difference of order unity and when this difference is introduced inside the generalized Ohm law.

Then, we propose a kinetic formulation for the Hall-MHD equations which contains as fluid closure different variants of the Hall-MHD model. The kinetic formulation consists of a Fokker-Planck equation for the ions and a set of fluid equations for the electrons coupled through quasineutrality. The Fokker-Planck operator models electron-ion collisions and contributes to relaxation of the velocities and the temperatures of both species to a common value. This kinetic model was in particular justified in [START_REF] Degond | Transport coefficients of plasmas and disparate mass binary gases[END_REF]. The MHD equations are obtained by taking the fluid moments of the ion Fokker-Planck equation and closing the resulting equations by a Maxwellian assumption. The resulting two-temperature resistive Hall MHD model consists of conservation equations for the density, momentum, energy and magnetic field combined with an evolution equation for the electron temperature and with the generalized Ohm law. The proposed hybrid ion-kinetic, electron-fluid model bear strong analogies with models used in the literature for numerical simulations such as [START_REF] Valentini | A hybrid-Vlasov model based on the current advance method for the simulation of collisionless magnetized plasma[END_REF], but the purpose is to highlight its mathematical structures. Indeed, a particularly interesting special case is when the electron and ion temperatures are equal to the same constant value (isothermal single temperature resistive Hall MHD). In this case, we can rephrase its kinetic formulation in the form of a coupled Fokker-Planck Faraday system, which exhibits an entropy dissipation identity. Surprisingly enough, the kinetic formulation of standard ideal MHD is deduced by neglecting the Hall term in Faraday's equation but keeping it in the kinetic equation.

The theoretical analysis focuses on the existence of global weak solutions for the incompressible viscous resistive Hall-MHD model written as follows:

∂ t u + u • ∇u + ∇p = (∇ × B) × B + ∆u, (1.1) ∇ • u = 0, (1.2) ∂ t B -∇ × (u × B) + ∇ × ((∇ × B) × B) = ∆B, (1.3) ∇ • B = 0, (1.4)
where u(x, t) and B(x, t) are the fluid velocity and magnetic field, depending on the spatial position x and the time t. The result is valid on a square domain Ω of R 3 with periodic boundary conditions. We stress that the important contribution of this work is the account of the last term of the left-hand side of (1.3), known as the Hall effect term.

The main theorem of this work is stated as follows:

Theorem 1.1 Let Ω = [0, 1] 3 . Assume that u 0 ∈ (L 2 (Ω)) 3 , B 0 ∈ (L 2 (Ω)) 3 with ∇•u 0 = 0, ∇ • B 0 = 0.
Then, there exists a global weak solution (u, B) for the Hall MHD problem (1.1), (1.4). Moreover, we have

(u, B) ∈ L ∞ ((0, T ), L 2 (Ω)) ∩ L 2 ((0, T ), H 1 (Ω)) and ∂ t u ∈ L 4 3 ((0, T ), H -1 (Ω)), ∂ t B ∈ L 4 3 ((0, T ), H -2 (Ω)).
Additionally, the following energy inequality holds:

d dt E(t) + ∇u 2 L 2 (Ω) + ∇B 2 L 2 (Ω) ≤ 0, (1.5) with E(t) = 1 2 u 2 L 2 (Ω) + B 2 L 2 (Ω) . (1.6)
In fact, the main difficulty of this work is concentrated in the treatment of the Hall term. So, we will show a preliminary result for the following Hall problem

∂ t B + ∇ × ((∇ × B) × B) = ∆B, (1.7) ∇ • B = 0, (1.8)
and will provide a detailed proof. The proof uses the particular structure of the Hall term which has zero contribution to the energy identity. The proof of the existence for the coupled system is then a direct consequence of the energy inequality (1.5) and will only be sketched. This result is up to our knowledge the first theoretical result for Hall MHD.

Global existence for standard viscous resistive incompressible MHD has been previously proved by Duvaut & Lions [START_REF] Duvaut | inéquations en thermoélasticité et magnétohydrodynamique[END_REF]. A bifurcation analysis of the Hall-MHD problem in view of the question of magnetic reconnection is performed in [START_REF] Homann | Bifurcation analysis of magnetic reconnection in Hall-MHD systems[END_REF]. Numerical methods for solving the Hall-MHD problems can be found e.g. in [START_REF] Arnold | A semi-implicit Hall-MHD solver using whistler wave preconditioning[END_REF][START_REF] Chacòn | A 2D high-β Hall MHD implicit nonlinear solver[END_REF][START_REF] Harned | Accurate semi-implicit treatment of the Hall effect in magnetohydrodynamic computations[END_REF][START_REF] Huba | Hall magnetohydrodynamics of reversed field current layers[END_REF][START_REF] Kazeminezhad | A discrete model for MHD incorporating the Hall term[END_REF] Remark 1.1 In a general domain Ω, the physically relevant boundary condition is the perfectly conducting wall boundary condition, which consists in assuming zero normal component of the B field and zero tangential component of the electric field. Here, the electric field is the quantity inside the curl operator, namely

E = (∇ × B) × B + ∇ × B.
This leads to the nonlinear boundary conditions

B • n = 0 on ∂Ω, n × (∇ × B) + n × ((∇ × B) × B) = 0 on ∂Ω.
From these conditions, we deduce that

B • (∇ × B) = 0,
which means that there is no helicity on the boundary. Because of the nonlinearity of this boundary condition, the methods developed below do not apply.

Finally, we discuss particular solutions of the Hall problem in the form of axisymmetric purely swirling magnetic fields and propose some regularization of the Hall equation. For axisymmetric purely swirling magnetic fields, the Hall problem reduces to a viscous Burger's equation. By neglecting the resistivity, the resulting inviscid Burger's equation shows shock wave solutions which are known in physics textbooks as KMC waves for Kingsep, Mokhov and Chukbar [START_REF] Kingsep | Nonlinear skin phenomenas in plasmas, Nonlinear and Turbulent Processes in Physics[END_REF] (see also [START_REF] Cassany | Analysis of the operating regimes of microsecond-conduction-time plasma opening switches[END_REF][START_REF] Dreher | Axisymmetric flows in Hall-MHD: a tendency towards finite-time singularity formation[END_REF]). They only exist if the Hall term is present. This effect also generates boundary layers which lead to nonlinear boundary conditions (see e.g. [START_REF] Méhats | A nonlinear oblique derivative boundary value problem for the heat equation. Part 1: basic results[END_REF]). Focusing on the non-resistive Hall problem itself, we propose a regularization consisting in restoring the displacement current in the Ampere equation. We then provide two equivalent formulations of this regularized problem which are obtained when either the current or the electric fields are eliminated from the system.

The organization of the paper is as follows. In Section 2, we propose a derivation of the Hall-MHD model from the two-fluids Euler-Maxwell model under suitable scaling hypotheses. Then, Section 3 is devoted to the presentation of the kinetic formulation of the Hall-MHD problem. Section 4 is focuses on the proof of the existence of global weak solutions for the incompressible viscous resistive Hall-MHD equations. Section 5 discusses the particular case of axisymmetric purely swirling magnetic fields and proposes a regularization of the Hall problem by means of a re-introduction of the displacement current in the Ampere equation. Finally, a conclusion is drawn in Section 6.

Derivation of the Hall-MHD equations from fluid equations

In this section, we briefly motivate the derivation of the model we are considering. For simplicity, we consider the compressible inviscid model and later on change to viscous incompressible flow. We start from the two-fluid isothermal Euler-Maxwell system for the electrons and ions, where we assume that the electron and ion temperatures are the same given constant:

∂ t n e + ∇ • (n e u e ) = 0, ( 2
∂ t n i + ∇ • (n i u i ) = 0, m i (∂ t (n i u i ) + ∇(n i u i ⊗ u i )) + ∇(n i T ) = = en i (E + u i × B) -e 2 ηn i n e (u i -u e ), c -2 ∂ t E -∇ × B = -µ 0 j, ǫ 0 ∇ • E = ρ, ∂ t B + ∇ × E = 0 , ∇ • B = 0, ρ = e(n i -n e ), j = e(n i u i -n e u e ).
where n e and n i are the electron and ion densities, u e and u i , the velocities, T , the common electron and ion temperature, m e and m i the masses. e denotes the elementary positive charge, and we assume singly charged positive ions. η is the resistivity due to the electronion collisions. E, B, ρ, j are respectively the electric field, the magnetic field, the charge density and the current density. ǫ 0 , µ 0 and c are respectively the vacuum permittivity, the vacuum permeability and the speed of light, related by the relation ǫ 0 µ 0 c 2 = 1. For simplicity, we assume monoatomic perfect gas equations of states for both the electrons and ions. We make the Boltzmann constant equal to unity which means that we measure temperatures in units of energy. The last terms at the right-hand sides of the second and fourth equations are the contributions of the electron-ion collisions to the momentum equation of each species. The two terms sum up to zero which expresses the conservation of total momentum in such collisions.

We introduce scaling units n 0 , u 0 , E 0 , B 0 , x 0 , t 0 , ρ 0 , j 0 for respectively the densities, velocities, electric field, magnetic field, space, time, charge and current. We assume that these units are related by the following relations:

x 0 = u 0 t 0 , u 0 = T m i , E 0 = u 0 B 0 , ρ 0 = en 0 .
The first relation means that we observe the system at the convection time scale. The second relation states that the convection velocity is that of the ion thermal speed. The third relation is typical of a MHD scaling and states that the main contribution to the electric field is induction due to the motion of the charged fluid. Finally, the last relation expresses the consistency between the density and charge units. Then, six dimensionless parameters appear:

ε 2 = m e m i , α 2 = eE 0 x 0 T , β = e 2 ηn 0 u 0 x 0 T γ = u 0 c , λ 2 = ε 0 T e 2 n 0 x 2 0 , η = j 0 en 0 u 0 .
which have the following interpretation. ε 2 is the electron to ion mass ratio and is very small. α 2 is the ratio of the electric energy to the thermal energy. β measures the relaxation frequency of the electron and ion velocities due to collisions. γ is the ratio of the fluid velocity to the speed of light. λ is the scaled Debye length and measures the closedness to quasi-neutrality. η is the ratio of the charge current scale to the electron or ion current scales. Since the charge current is the difference of these two particle currents, it may be much smaller than any of them due to charge neutrality. Therefore, the scale ratio η may be either O(1) or ≪ 1 according to the situations. The dimensionless two-fluids Euler-Maxwell system is written:

∂ t n e + ∇ • (n e u e ) = 0, ε 2 (∂ t (n e u e ) + ∇(n e u e ⊗ u e )) + ∇(n e T ) = = -α 2 n e (E + u e × B) -βn e n i (u e -u i ), ∂ t n i + ∇ • (n i u i ) = 0, ∂ t (n i u i ) + ∇(n i u i ⊗ u i ) + ∇(n i T ) = = α 2 n i (E + u i × B) -βn i n e (u i -u e ), γ 2 ∂ t E -∇ × B = - γ 2 η α 2 λ 2 j, α 2 λ 2 ∇ • E = ρ, ∂ t B + ∇ × E = 0 , ∇ • B = 0, ρ = n i -n e , j = 1 η (n i u i -n e u e ).
The compressible MHD equation corresponds to the simultaneous independent four limits 1. ε 2 → 0: this corresponds to the neglect of the convection term in the electron momentum equation. The resulting equation is usually referred to as the generalized Ohm's law.

2. λ 2 → 0. This gives rise to quasineutrality, i.e. the fact that the local electron and ion densities are everywhere the same. We now denote by n their common value:

n e = n i = n.
3. γ 2 → 0 while keeping γ 2 η α 2 λ 2 = 1. This leads to the neglect of the displacement current in Ampere's equation and gives rise to the standard magnetostatic Ampere law.

The resulting system is the so-called compressible isothermal resistive Hall-MHD equations. Denoting by u the ion velocity, this system is written

∂ t n + ∇ • (nu) = 0, ∂ t (nu) + ∇(nu ⊗ u) + ∇(2nT ) = α 2 η j × B, (2.3) ∇ × B = j, ∂ t B + ∇ × E = 0, ∇ • B = 0, j = 1 η n(u -u e ), (2.4)
E + u × B = - T α 2 ∇(ln n) + η j × B n + βη α 2 j, (2.5)
where we highlight the momentum conservation eq. ( 2.3), the current equation (2.4) and generalized Ohm's law (2.5). Note that the T ∇ln n term at the right-hand side of (2.5) has no contribution since the curl operator in the Faraday equation cancels it. However, this cancellation is no more true in the general gas dynamics case because ∇pe n may not be a gradient in general, where p e is the electron pressure.

In all what follows, we assume α 2 η = 1 in order to keep the Lorentz force term in (2.3) of order 1 in all the various scalings below. Then, eqs (2.3), (2.4) and (2.5) are written :

∂ t (nu) + ∇(nu ⊗ u) + ∇(2nT ) = j × B, (2.6) 1 α 2 j = n(u -u e ), (2.7)
E + u × B = 1 α 2 -T ∇(ln n) + j × B n + β α 4 j, (2.8)
the other equations being unchanged. There are only two dimensionless parameters left:

1 α 2 and β α 4 and they only appear in (2.7) and in (2.8). So, the various types of MHD model correspond to the various choices of scalings for these two parameters. In particular, we have 1. If both 1 α 2 → 0 and β α 4 → 0, then the generalized Ohm's law reduces to the standard ideal Ohm's law while the electron and ion velocities become identical:

E + u × B = 0, u e = u.
This yields ideal MHD.

If 1

α 2 → 0 but β α 4 → 1, then the resistive term in the generalized Ohm's law is kept but the electron and ion velocities are still identical:

E + u × B = j, u e = u.
This gives rise to resistive MHD.

If 1

α 2 → 1 but β α 4 → 0, then, the ion and electron velocities differ and additionally, the generalized Ohm's law has the form:

E + u × B = -T ∇(ln n) + j × B n .
As already mentioned, the first term at the right-hand side has no contribution. The second one is the Hall term. This gives rise to the Hall MHD.

4. Finally, if both 1 α 2 → 1, β α 4 → 1, then, the ion and electron velocities differ and both the Hall and resistive terms appear.

E + u × B = -T ∇(lnn) + j × B n + j.
Our study takes place in the context of the last regime, where both the resistive and Hall terms are equally important. Additionaly, we assume incompressible viscous fluid motion. In this case, the Hall MHD system can be written according to (1.1)- (1.4).

We note that it is easy to extend this system to the viscous isentropic resistive compressible Hall MHD as follows (assuming all the physical constants equal to 1 except the viscosity here denoted by ν):

∂ t n + ∇ • (nu) = 0, (2.9) ∂ t (nu) + ∇(nu ⊗ u) + ∇p(n) = (∇ × B) × B + ν∇ • (∇u + (∇u) T ),
(2.10)

∂ t B + ∇ × B × u + (∇ × B) × B n = -∇ × (∇ × B), (2.11)
For this system, we have the following magneto-helicity conservation relation:

d dt R 3 B • A dx + 2 R 3 B • (∇ × B) dx = 0, (2.12)
where A such that B = ∇ × A is any vector potential of B. To prove this relation, we note, using Green's formula, that

d dt R 3 B • A dx = 2 R 3 A • B t dx.
Then, taking the scalar product of (2.11) with 2A easily gives the result.

Remark 2.1

The viscosity term at the right-hand side of (2.10) involves the rate of strain tensor σ(u) = ∇u+(∇u) T . In usual gas dynamics, the viscosity term involves the traceless rate of strain tensor σ 0 (u) = ∇u+(∇u) T -(2/3)(∇•u)Id. However, a careful computation of the viscosity in the case of isothermal gas dynamics shows that, in this case, the right tensor is the full rate-of-strain tensor σ(u) and not its trace-free counterpart σ 0 (u).

Derivation of the Hall-MHD equations from kinetic equations

In this section, we provide a kinetic formulation of the Hall MHD problem. We start from a kinetic equation for the ion distribution function f (x, v, t) of the plasma, where x is the position, v the velocity and t the time. This distribution function is a solution of the following kinetic equation

∂ t f + v • ∇ x f + e m (E + v × B) • ∇ v f = Q(f ) , (3.1)
where e is the positive ion charge, supposed equal to the absolute value of the elementary charge, m is their mass, E(x, t) and B(x, t) are the electric and magnetic fields, and Q(f ) is the collision operator for electron-ion collisions. We respectively introduce the ion density n, mean velocity u and energy W by

n = f dv , nu = f v dv , W = f m |v| 2 2 dv,
and the temperature by

3 2 nT = W - 1 2 n|u| 2 = f m |v -u| 2 2 dv.
We assume that the electrons are described by their fluid quantities, namely their density n e (x, t), their fluid velocity u e (x, t) and their temperature T e (x, t). The use of a fluid model for the electrons while the ions are treated kinetically can be justified by the small electron to ion mass ratio. A formally rigorous justification of this can be found e.g. in [START_REF] Degond | Transport coefficients of plasmas and disparate mass binary gases[END_REF]. The electron density is supposed equal to the ion density by quasineutrality:

n e = n . (3.2)
The electron momentum conservation equation, when the transport term is neglected due to their small mass, gives rise to the generalized Ohm's law (see (2.2) where all terms in factor of m e are set to zero):

∇ x (nT e ) + en(E + u e × B) = eηnj , (3.3)
where m e is the electron mass, j is the current density and η is the resistivity. The current density is given by:

j = en(u -u e ) . (3.4)
Since the electron mass is neglected and assuming monoatomic gas equation of state, the electron energy W e can be expressed in terms of the electron temperature T e by

W e = 3 2 nT e .
For the same reason, the mass is neglected in the electron energy flux, which reads 5 2 nT e u e . For simplicity, we consider a model electron-ion collision operator as follows:

Q(f ) = e 2 ηn m ∇ v • ((v -u e )f + T e m ∇ v f ) . (3.5)
The first term expresses the relaxation of the ion velocity to the electron one, while the second one expresses the relaxation of the ion temperature to the electron one. More realistic expressions of the electron-ion collision operator can be found in the literature (see e.g. [START_REF] Braginskii | Transport processes in a plasma[END_REF]), but this model is chosen for the sake of simplicity of exposition. The rate of change of the ion momentum is given by:

Q(f ) mv dv = -eηnj,
and is the opposite of the right-hand side of (3.3), which is consistent with the total momentum conservation of the electron-ion collisions. The rate of change of the ion energy is given by

Q(f ) m|v| 2 2 dv = -eηnj • u + 3 e 2 ηn 2 m (T e -T ).
By total energy conservation in electron-ion collisions, the rate of change of the electron total energy is the opposite. Then, the electron energy conservation equation reads: In this equation, the first term of the right-hand side is the work done by the electrons in the Lorentz force, while the last two terms are due to the electron-ion collisions. Taking the scalar product of (3.3) by u e and subtracting it to (3.6) leads to

∂ t 3 2 nT e + ∇ x • 5
(∂ t + u e • ∇ x ) 3 2 nT e + 5 2 nT e ∇ x • u e = η|j| 2 + 3 e 2 ηn 2 m (T -T e ). (3.7)
The first term of the right-hand side is Joule heating of the electrons, while the second term is the electron temperature relaxation to the ion temperature.

The magnetic field evolves according to the Faraday equation

∂ t B + ∇ x × E = 0 , (3.8)
and the current is linked to the magnetic field by Ampere's law:

∇ x × B = µ 0 j , (3.9)
As in the previous section, the displacement current has been neglected.

As a summary, the considered kinetic model is as follows:

∂ t f + v • ∇ x f + e m (E + v × B) • ∇ v f = e 2 ηn m ∇ v • ((v -u e )f + T e m ∇ v f ) , (3.10)
∇ x (nT e ) + en(E + u e × B) = eηnj , (3.11)
(∂ t + u e • ∇ x ) 3 2 nT e + 5 2 nT e ∇ x • u e = η|j| 2 + 3 e 2 ηn 2 m (T -T e ) , (3.12
)

∂ t B + ∇ x × E = 0 . (3.13) ∇ x × B = µ 0 j , (3.14) j = en(u -u e ) , (3.15
)

n = f dv , nu = f v dv , 3 2 nT = f m |v -u| 2 2 dv. (3.16)
Now, we link this system to Hall-MHD by taking the moments of the ion kinetic equation. Integrating (3.10) with respect to v after premultiplying it successively by 1, mv or m |v| 2 2 , we get the following ion mass, momentum and energy balance equations:

∂ t n + ∇ x • (nu) = 0 , (3.17) m (∂ t (nu) + ∇ x • (nu ⊗ u)) + ∇ x • P = en(E + u × B) -eηnj , (3.18)
∂ t W + ∇ x • (W u + Pu + q) = enE • u -eηnj • u + 3 e 2 ηn 2 m (T e -T ) , (3.19)
where P and q are the stress tensor and heat flux vector, given by:

P = m f (v -u) ⊗ (v -u) dv, q = m 2 f (v -u)|v -u| 2 dv . ( 3
(∇ x × B) × B = ∇ x • (B ⊗ B) -∇ x (|B| 2 /2
), the total fluid momentum balance can be written in conservative form:

∂ ∂t (mnu) + ∇ x • (mnu ⊗ u - 1 µ 0 B ⊗ B + P tot ) = 0 , (3.22)
where the total pressure tensor P tot is written

P tot = P + (nT e + |B| 2 2µ
f = W + W e : ∂ ∂t W f + ∇ x • ((W e + p e )u e + W u + Pu + q) = E • j . (3.24)
The Faraday equation (3.13) implies that ∂ ∂t

|B| 2 2µ 0 + 1 µ 0 B • (∇ x × E) = 0 . (3.25)
By adding (3.24) and (3.25) and using Ampere's law (3.14), the total energy (which is the sum of the total fluid energy and the magnetic energy) W tot = W f + |B| 2 2µ 0 satisfies the following conservation law:

∂ ∂t W tot + ∇ x • ((W e + p e )u e + W u + Pu + q + 1 µ 0 E × B) = 0 , (3.26)
The MHD equations are obtained through the closure assumptions that P = nT Id, q = 0, which can be justified e.g. by a Maxwellian closure, i.e. assuming that f = M n,u,T with

M n,u,T = n 2πT m 3/2 exp - m|v -u| 2 2T .
The Maxwellian closure itself can be justified if ion-ion collisions are strong enough to relax the distribution f quickly to M n,u,T . However, in many instances, the Maxwellian closure is used in spite of not being fully justified. In this case, we obtain the 2-temperature, Hall, resistive compressible MHD equations which are as follows:

∂ t n + ∇ x • (nu) = 0 , (3.27) ∂ t (mnu) + ∇ x • (mnu ⊗ u - 1 µ 0 B ⊗ B + (n(T e + T ) + |B| 2 2µ 0 )Id) = 0 , (3.28) ∂ t W + 3 2 nT e + |B| 2 2µ 0 + ∇ x • (W u + nT u + 5 2 nT e u e + 1 µ 0 E × B) = 0 , (3.29) ∂ t B + ∇ x × E = 0 , (3.30) ∇ x (nT e ) + en(E + u e × B) = eηnj , (3.31) (∂ t + u e • ∇ x ) 3 2 nT e + 5 2 nT e ∇ x • u e = η|j| 2 + 3 e 2 ηn 2 m (T -T e ) , (3.32)
∇ x × B = µ 0 j , (3.33) j = en(u -u e ).
(3.34)

The first four equations are the basic conservation laws of mass, total momentum, total energy and magnetic field. Eq. (3.31) is the generalized Ohm law and provides the expression for E. Eq. ( 3.32) provides the evolution equation for T e . Finally, Ampere's eq. (3.33) defines j and (3.34) defines u e . The fact that u e = u gives rise to the Hall term. Now, if T = T e , then, eq. ( 3.32) can be removed and one gets the single temperature Hall resistive compressible MHD equations. Another simplification is to suppose that the common ion and electron temperatures are constant (isothermal assumption). In this case, the total energy equation (3.29) is a consequence of the Faraday and momentum eqs. (3.30), (3.28) and can be removed. Then, we find the model of section 2.

A kinetic formulation of the model considered in section 2 is also obtained from the kinetic model (3.10)- (3.16) by supposing that the electron and ion temperatures are the same constant T . This kinetic model is written below:

∂ t f + v • ∇ x f + e m (E + v × B) • ∇ v f = e 2 ηn m ∇ v • ((v -u e )f + T m ∇ v f ) , (3.35) T ∇ x n + en(E + u e × B) = eηnj , (3.36) ∂ t B + ∇ x × E = 0 , (3.37) ∇ x × B = µ 0 j , (3.38) j = en(u -u e ) .
(3.39)

By eliminating E, j and u e respectively using (3.36), (3.38), (3.39), and after some easy algebra, we find the following coupled Fokker-Planck Faraday system for f and B:

∂ t f + v • ∇ x f + e m (v -u) × B + 1 µ 0 en (∇ x × B) × B - T e ∇ x ln n • ∇ v f = = e 2 ηn m ∇ v • ((v -u)f + T m ∇ v f ) , (3.40) ∂ t B + ∇ x × B × u + 1 µ 0 en (∇ x × B) × B = - η µ 0 ∇ x × (∇ x × B) . (3.41)
For the kinetic eq. (3.40), we have the following entropy dissipation identity associated to the entropy H(f ):

dH(f ) dt + e 2 η T (x,v)∈R 3 ×R 3 n (v -u)f + T m ∇ v f 2 f dx dv = 0, H(f ) = (x,v)∈R 3 ×R 3 f ln f dx dv.
Particle-in-Cell simulations of this model can be found in [START_REF] Degond | Simulation of non-equilibrium plasmas with a numerical noise-reduced particle-in-cell method[END_REF]. Now, taking the first two moments of (3.40), we get the isothermal resistive Hall MHD model which was written in dimensionless form at (2.9)-(2.11) (taking ν = 0 and p(n) = 2T n). It is interesting to note that in order to get a kinetic model for the standard resistive isothermal MHD equations, we need to neglect the Hall term 1 µ 0 en (∇ x × B) × B in the Faraday eq. (3.41), but we must retain the corresponding term in the kinetic equation (3.40). Therefore, a kinetic formulation for standard resistive isothermal MHD is as follows:

∂ t f + v • ∇ x f + e m (v -u) × B + 1 µ 0 en (∇ x × B) × B - T e ∇ x ln n • ∇ v f = = e 2 ηn m ∇ v • ((v -u)f + T m ∇ v f ) , (3.42)
∂ t B + ∇ x × (B × u) = - η µ 0 ∇ x × (∇ x × B) . (3.43)
We can also get a kinetic formulation for the ideal isothermal MHD equations by neglecting the resistive term η µ 0 ∇ x × (∇ x × B) at the right-hand side of the Faraday eq. (3.43). The collision term at the right-hand side of (3.42) has no contribution in the mass and momentum conservation equations and can be either kept or removed without modifying the corresponding balance equations.

4 Existence result for the incompressible, viscous, resistive Hall-MHD

In the paper, we use the function spaces

H per (Ω) = {B ∈ (H 1 (Ω)) 3 | ∇ • B = 0 on Ω with periodicity conditions}.
We use A, B = Ω A • B dx, for any pair A, B ∈ (L 2 (Ω)) 3 . We remark that

∇B (L 2 (Ω)) 9 = ∇ × B (L 2 (Ω)) 3 . (4.1)
To prove theorem 1.1, we just focus on the Hall problem itself. The extension to the full Hall-MHD problem is explained at the end of the section. So, our goal is now to prove the following existence theorem for the Hall problem. We introduce the following weak formulations of the Hall MHD problem:

Find B ∈ L ∞ (0, T, L 2 (Ω)) ∩ L 2 ((0, T ), H per (Ω)) such that for any A ∈ H per (Ω): A, ∂ t B + ∇ × A, (∇ × B) × B + ∇ × A, ∇ × B = 0. (4.2) Theorem 4.1 Assume that B 0 ∈ (L 2 (Ω)) 3 with ∇ • B 0 = 0.
Then, there exists a global weak solution B for the Hall problem (1.7), (1.8). Moreover, we have

B ∈ L ∞ ((0, T ), L 2 (Ω)) ∩ L 2 ((0, T ), H per (Ω)) and ∂ t B ∈ L 4 3 ((0, T ), H -2 (Ω)).
The proof is based on the construction of an approximate solution by Galerkin's method. Uniform a priori bounds on these approximate solutions will allow us to pass to the limit thanks to standard compactness arguments.

We define the Fourier basis (φ k ) k∈Z 3 , with φ k = e 2πik•x . We denote by |k| = |k 1 |+|k 2 |+ |k 3 |. The Fourier basis provides a complete ortho-normal basis of L 2 (Ω). We denote by P N the projection onto the Fourier basis with index |i| ≤ N. We define an approximate solution B N : (0, T ) → H per (Ω) of problem (1.1), (1.3) of the form:

B N (t) = i∈Z 3 , |i|≤N B N i (t)ϕ i (4.3)
with the divergence free constraint i • B N i (t) = 0, satisfying

A N , ∂ t B N (t) + ∇ × A N , (∇ × B N (t)) × B N (t) + ∇ × A N , ∇ × B N (t) = 0, (4.4)
for all A N of the form

A N = N i=1 A N i ϕ i , i • A N i = 0,
and with initial condition B N (0) = P N B 0 .

Lemma 4.2 There exists a global in time solution B N (t) of (4.4) which is uniformly

(independently of N) bounded in L ∞ ((0, T ), L 2 (Ω)) ∩ L 2 ((0, T ), H 1 (Ω)) such that ∂ t B N is uniformly bounded in L 4 3 ((0, T ), H -2 (Ω)).
Proof. Throughout of the proof, C denotes a generic constant. We take A = B N in (4.4) and get:

B N , ∂ t B N + ∇ × B N 2 (L 2 (Ω)) 3 = 0. (4.5)
Thus, we have:

B N (t) 2 L 2 (Ω) + 2 t 0 ∇B N (s) 2 (L 2 (Ω)) 9 ds = B N 0 2 L 2 (Ω) . (4.6)
This shows the uniform bound of B N in L ∞ ((0, T ), L 2 (Ω)) ∩ L 2 ((0, T ), H 1 (Ω)).

To show the time regularity, we take A ∈ H 2 (Ω) 3 , such that ∇ • A = 0 and take

A N = P N A. A, ∂ t B N = A N , ∂ t B N = -∇ × A N , (∇ × B N ) × B N -∇ × A N , ∇ × B N . (4.7)
We recall that

(∇ × B) × B = ∇ • (B ⊗ B) -∇ |B| 2 2 .
(4.8)

Introducing (4.8) into (4.7) and noting that, by Green's formula, a gradient and a curl are orthogonal in L 2 , we find

A, ∂ t B N = ∇ × A N , ∇ • (B N ⊗ B N ) + ∇ × A N , ∇ × B N . (4.9)
To estimate the first term at the right-hand side of (4.9), we use Green's formula and get:

∇ × A N , ∇ • (B N ⊗ B N ) = -B N , (B N • ∇)∇ × A N .
Using Hölder's inequality, and remarking that A N H 2 (Ω) ≤ A H 2 (Ω) , we get

| B N , (B N • ∇)∇ × A N | ≤ C B N 2 L 4 (Ω) A N H 2 (Ω) ≤ C B N 2 L 4 (Ω) A H 2 (Ω) . (4.10)
By Gagliardo-Nirenberg's inequality [START_REF] Evans | Partial differential equations[END_REF],

B N 2 L 4 (Ω) ≤ C ∇B N 3 2 L 2 (Ω) ( B N 1 2 L 2 (Ω) + 1),
and using the uniform L 2 bound, we can write:

| B N , (B N • ∇)∇ × A N | ≤ C ∇B N 3 2 L 2 (Ω) A H 2 (Ω)
. The second term is simply estimated by using Cauchy-Schwartz inequality:

| ∇ × A N , ∇ × B N | ≤ C ∇B N L 2 (Ω) A N H 1 (Ω) ≤ C ∇B N L 2 (Ω) A H 1 (Ω)
Collecting these estimates, we obtain:

| A, ∂ t B N | ≤ C ∇B N 3 2 L 2 (Ω) A H 2 (Ω) + ∇B N L 2 (Ω) A H 1 (Ω) ≤ C ∇B N 3 2 L 2 (Ω) + 1 A H 2 (Ω) .
Therefore,

∂ t B N H -2 (Ω) ≤ C ∇B N 3 2 L 2 (Ω) + 1 . Thus: ∂ t B N 4 3 H -2 (Ω) ≤ C ∇B N 2 L 2 (Ω) + 1 .
Thanks to (4.6), the right-hand side is integrable on (0, T ) and we get

T 0 ∂ t B N 4 3 H -2 (Ω) dt ≤ C, (4.11)
which ends the proof of Lemma 4.2.

Thanks to this Lemma, we can proceed to the Proof of Theorem 4.1. According to Lemma 4.2, the sequence

(B N ) N ∈N is uniformly bounded in L 2 ((0, T ), H 1 (Ω)) with (∂ t B N ) N ∈N uniformly bounded in L 4
3 ((0, T ), H -2 (Ω)). Consequently, by virtue of Lions-Aubin Lemma [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF], (B N ) N ∈N is compact in L 2 ((0, T ), L 2 (Ω)). Therefore, there exists a subsequence (B N k ) k∈N and a function B in L 2 ((0, T ), H 1 (Ω)) with ∂ t B in L 4 3 ((0, T ), H -2 (Ω)), such that:

           B N k ⇀ B weak star in L ∞ ((0, T ), L 2 (Ω)) B N k ⇀ B weakly in L 2 ((0, T ), H 1 (Ω)) B N k → B strongly in L 2 ((0, T ), L 2 (Ω)) ∂ t B N k ⇀ ∂ t B weakly in L 4 3 ((0, T ), H -2 (Ω)) (4.12)
We take A ∈ H 3 (Ω) 3 , such that ∇ • A = 0, take A N k = P N k A in (4.4) and integrate it with respect to time. We have: 3 . Now, thanks to the convergences (4.12), we can take the limit N k → ∞ and get that B satisfies

A N k , B N k (t) -A N k , B N k (0) + t 0 ∇ × A N k , (∇ × B N k ) × B N k ds + + t 0 ∇ × A N k , ∇ × B N k ds = 0. Thanks to the Sobolev imbedding in dimension 3, ∇ × A ∈ L ∞ (Ω) 3 and ∇ × A N k → ∇ × A strongly in L ∞ (Ω)
A, B(t) -A, B(0) + t 0 ∇ × A, (∇ × B) × B ds + + t 0 ∇ × A, ∇ × B ds = 0,
which is a weak solution of the Hall problem.

Proof of theorem 1.1. We apply the same Galerkin construction (u N , B N ) for the coupled system (1.1), (1.3) as we did for the Hall system and use the energy identity for the Galerkin approximation:

u N (t) 2 L 2 (Ω) + B N (t) 2 L 2 (Ω) + 2 t 0 ( ∇u N (s) 2 (L 2 (Ω)) 9 + ∇B N (s) 2 (L 2 (Ω)) 9 ) ds = = u N 0 2 L 2 (Ω) + B N 0 2 L 2 (Ω) .
Then the same proof can be reproduced. We just note that the time regularity of u can be improved because there is no Hall term in the velocity equation. Therefore, we find

u t ∈ L 4 3 ((0, T ), H -1 (Ω)) while B t ∈ L 4 3 ((0, T ), H -2 (Ω)).

Axisymmetric flows, KMC waves and Maxwell regularization of the non-resistive Hall problem

In this section, we assume axisymmetric B field. Let x = (x, y, z) a coordinate system, where x is the symmetry axis. Axisymmetry about the x-axis means that, given any rotation R about this axis, the field B satisfies: B(Rx) = RB(x). We use (x, r, θ) the cylindrical coordinates of x and (e x , e r , e θ ) as associated local basis. Then, using a representation of axisymmetric, divergence-free fields given in [START_REF] Liu | Characterization and regularity for axisymmetric solenoidal vector fields with application to Navier-Stokes equation[END_REF], we can write

B = be θ + ∇ × (ψe θ ),
where the scalar functions b and ψ are functions of (x, r). We note the simple formulas [START_REF] Liu | Characterization and regularity for axisymmetric solenoidal vector fields with application to Navier-Stokes equation[END_REF]:

∇ × B = -Lψe θ + ∇ × (be θ ), L = ∇ 2 (x,r) - 1 r 2 , ∇ × (∇ × B) = -Lb + ∇ × (je θ ), j = -Lψ, ∇ × ((∇ × B) × B) = j r , rψ x,r - b r , rb x,r e θ + +∇ × 1 r 2 {rb, rψ} x,r e θ ,
where the Poisson Bracket {a, b}

x,r = ∂ x a ∂ r b -∂ x b ∂ r a.
In these coordinates, the Hall problem (1.7), (1.8) is written:

ψ t + 1 r 2 {rb, rψ} = Lψ, b t + j r , rψ - b r , rb = Lb.
If initially ψ 0 = 0 (which means that the B field is purely swirling), then ψ ≡ 0 at all times and the b equation reduces to the simple viscous Burger's equation

b t - 2 r b b x = Lb.
In the inviscid case, the Burger's equation has solutions in the form of propagating shock waves. In the context of Hall MHD, these waves are known as the KMC waves, for Kingsep, Mokhov and Chukbar [START_REF] Kingsep | Nonlinear skin phenomenas in plasmas, Nonlinear and Turbulent Processes in Physics[END_REF] (see also [START_REF] Cassany | Analysis of the operating regimes of microsecond-conduction-time plasma opening switches[END_REF]). They only exist if the Hall term is present.

Here, we focus on the non-resistive Hall problem and formulate it as a limit of a so-called Maxwell regularization. The non-resistive Hall problem is written

∂ t B + ∇ × E = 0, ∇ × B = j, E = j × B.
In the case of axisymmetric purely swirling B-fields, this yields the inviscid Burger's equation as seen above. We now consider a regularization of this problem by restoring the displacement current in the Ampere equation. This yields the problem

∂ t B + ∇ × E = 0,
(5.1) -ε∂ t E + ∇ × B = j, (5.

2) E = j × B, (5.3) where ε ≪ 1 is the Maxwell regularization parameter. We now investigate two formulations of the Maxwell-regularized problem.

1. The (B, j) formulation. This formulation consists in classically eliminating E for j and gives ∂ t B + ∇ × (j × B) = 0, -ε∂ t (j × B) + ∇ × B = j.

When ε → 0, it clearly tends to the non-resistive Hall problem.

2. The (B, E) formulation. This formulation consists in eliminating j for E. From (5.3), we can write

j = B × E |B| 2 + λB,
where λ is a priori unknown. However, we also have the geometric constraint

E • B = 0, (5.4)
and λ can be viewed as a Lagrange multiplier of this constraint. Indeed, from

-ε∂ t E + ∇ × B = B × E |B| 2 + λB.
and taking the scalar product of this equation with B, we get

-ε(∂ t E) • B + (∇ × B) • B = λ|B| 2 .
But, differentiating the constraint (5.4) gives

(∂ t E) • B = -(∂ t B) • E = (∇ × E) • E.
Therefore,

λ = 1 |B| 2 (-ε(∇ × E) • E + (∇ × B) • B).
Finally, the Maxwell-regularized Hall problem in the (E, B) formulation reads

∂ t B + ∇ × E = 0, -ε∂ t E + ∇ × B = 1 |B| 2 (B × E + (-ε(∇ × E) • E + (∇ × B) • B) B).
In this system, the constraint (5.4) is satisfied as soon as it is satisfied at t = 0. The fact that the limit ε → 0 of this problem leads to the Hall problem is no more so obvious. Indeed, the limit ε → 0 leads to the following problem: and recover the Hall problem as the formal limit of the Maxwell-regularized system. It is an interesting problem to determine if this limit can be made rigorous.

∂ t B + ∇ × E = 0, ∇ × B = 1 |B| 2 (B × E + ((∇ × B) • B) B), ∂ t (E • B) = 0.

Conclusion

In this paper, we have derived and analyzed the Hall-MHD model. First, the model has been derived from a scaling limit of a two-fluids Euler-Maxwell system, under suitable scaling assumptions. Then, a derivation of the Hall-MHD from a kinetic model consisting of a Fokker-Planck equation for the ions coupled with fluid electrons has been realized. In the analysis section, we have proved the existence of global weak solutions for the incompressible viscous resistive Hall-MHD problem. The proof relies strongly on the skewsymmetric structure of the Hall term, which does not affect the energy estimate. This work shows that maintaining this structure is crucial for the well-posedness of the problem and is likely to be crucial as well for the derivation of stable numerical approximations. The question of the perfectly conducting boundary condition will be also investigated in future work. In the last section, we have reviewed some aspects of axisymmetric, purely swirling magnetic fields and proposed a regularization of the Hall problem by reintroducing the displacement current in Ampere's equation. Future work will be devoted to the analysis of the axisymmetric case and to the Maxwell regularization of the Hall problem, as well as to the investigation of the kinetic formulations of the Hall-MHD problem.

2

 2 nT e u e = -enu e • E + eηnj • u + 3 e 2 ηn 2 m (T -T e ). (3.6)

. 20 )

 20 By combining(3.18) with the generalized Ohm law(3.11), we obtain the total fluid momentum balancem ∂ ∂t (nu) + ∇ x • (nu ⊗ u) + ∇ x • (P + nT e Id) = j × B .(3.21)Using Ampere's equation(3.14), the fact that ∇ x • B = 0 and the vector identity

  The second equation is equivalent toB × [(∇ × B) × B -E] = 0.If (E • B)| t=0 = 0, then E • B ≡ 0 for all times; Then, we can invert this equation intoE = (∇ × B) × B,

  .1) m e (∂ t (n e u e ) + ∇(n e u e ⊗ u e )) + ∇(n e T ) = = -en e (E + u e × B) -e 2 ηn e n i (u e -u i ), (2.2)

Acknowledgements: This work has been supported by the 'Fondation Sciences et Technologies pour l'Aéronautique et l'Espace', in the frame of the project 'Plasmax' (contract # RTRA-STAE/2007/PF/002) who supported the fourth author's stay in Toulouse. This work has also been supported by the 'Commissariat à l' Énergie Atomique (CEA)' in the frame of the contract 'PICCADI' (contract # 4600192523). The second and fourth authors wish to acknowledge the hospitality of Mathematics Department and Mathematical

where this research was completed. The research of J.-G. L. was partially supported by NSF grant DMS 10-11738. The second author wishes to thank F. Deluzet, G. Gallice and C. Tessieras for fruitful discussions.


Julien Moras

Véronique Cherfaoui

Philippe Bonnifait

Credibilist Occupancy Grids for Vehicle Perception in Dynamic Environments

In urban environments, moving obstacles detection and free space determination are key issues for driving assistance systems and autonomous vehicles. When using lidar sensors scanning in front of the vehicle, uncertainty arises from ignorance and errors. Ignorance is due to the perception of new areas and errors come from imprecise pose estimation and noisy measurements. Complexity is also increased when the lidar provides multi-echo and multi-layer information. This paper presents an occupancy grid framework that has been designed to manage these different sources of uncertainty. A way to address this problem is to use grids projected onto the road surface in global and local frames. The global one generates the mapping and the local one is used to deal with moving objects. A credibilist approach is used to model the sensor information and to do a global fusion with the worldfixed map. Outdoor experimental results carried out with a precise positioning system show that such a perception strategy increases significantly the performance compared to a standard approach.

I. INTRODUCTION

There is a strong interest in the field of intelligent systems and robotics to develop autonomous vehicles with the ability to navigate in challenging environments. Indeed, several large-scale experiments like the DARPA Urban Challenge [START_REF] Buehler | The DARPA Urban Challenge: Autonomous Vehicles in City Traffic, George Air Force Base[END_REF] or the new challenge VIAC [START_REF] Broggi | Sensing requirements for a 13,000 km intercontinental autonomous drive[END_REF] (4 vehicles from Italy to China) have shown impressive results. However, autonomous driving in urban areas remains a problem and needs scene understanding capabilities to predict the short-term evolution for reliable navigation. For road navigation applications, safety is one of the fundamental requirements. A mere detection of the obstacles is not enough: a fine characterization of the navigable space is necessary to plan a safe motion of the robot [START_REF] Himmelsbach | Autonomous off-road navigation for mucar-3[END_REF]. For a car-like robot, navigable space means all the drivable space in the surrounding of the vehicle and obstacle means a insurmountable part of the space. This information has to be provided in real time with high level of confidence.

Some sensors providing high density information, for example the Velodyne [START_REF] Himmelsbach | Fast segmentation of 3d point clouds for ground vehicles[END_REF], are able to attempt the necessary level of confidence for navigation. But it requires perception systems that are capable of handling its high bandwidth or processing requirements. However, an other solution consists in using sensor providing sparse data. Representation of the environment is obtained exploiting temporal redundancy and temporal complementary of data.

In robotics, mapping (M) is coupled generally with localization in the Simultaneous Localization And Mapping (SLAM) problem [START_REF] Zhao | Slam in a dynamic large outdoor environment using a laser scanner[END_REF], [START_REF] Thrun | Probabilistic Robotics (Intelligent Robotics and Autonomous Agents)[END_REF]. Two main approaches exist: the The authors are with Heudiasyc,CNRS UMR 6599, Université de Technologie de Compiègne, 60200 Compiègne, France Contact author email: julien.moras@hds.utc.fr feature-based approach which tries to map with a predetermined set of shapes (segment, arc, etc...) [START_REF] Jeong | Visual slam with line and corner features[END_REF] and the grid approach which is based on a discrete space representation [START_REF] Garcia-Favrot | Laser scanner based slam in real road and traffic environment[END_REF]. SLAM methods are sensitive to the presence of moving objects in the scene. Indeed, algorithms are based on the temporal consistency of the map. A way to address this issue is to perform detection and Mobile Objects Tracking (MOT) [START_REF] Wang | Simultaneous localization, mapping and moving object tracking[END_REF]. Usually, a detection and tracking mechanism is used in the feature-based framework. It works in 3 steps: clustering of raw data for object detection, data association and temporal fusion [START_REF] Bar-Shalom | Multitarget-Multisensor tracking : Applications and Andvances[END_REF]. The main sources of error are in the clustering phase and in the association step [START_REF] Blackman | Design and Analysis of Modern Tracking Systems[END_REF]. Some recent works use parametric models to follow the variation of the shape of the objects [START_REF] Petrovskaya | Model based vehicle detection and tracking for autonomous urban driving[END_REF], [START_REF] Fayad | Tracking objects using a laser scanner in driving situation based on modeling target shape[END_REF].

This paper presents an approach dedicated to mapping and detection of moving obstacles using a sensor observing just a part of the environment. This kind of range sensor provides two types of information: the presence of echoes in the lineof-sight, and their relative position. The main contributions of this paper are the definition of a new sensor model adapted to a grid representation and a method to extract moving objects while performing the mapping using a credibilist fusion. As a clustering step is no longer necessary, there is no assumption on the shape of the objects, making it suitable for a wide range of urban objects (pedestrians, vehicles, bicycles ,...). Moreover, partially observed objects because of occlusions or sensor field of view limitation can be more precisely localized. The accurate pose of the robot is not a study of this paper and is provided by a specific module.

The paper is organized as follows. In section II, we present the occupancy grid. We then introduce the credibilist grids proposed in this work. Afterward, the sensor model is presented in section III, followed by the fusion processing stage in section IV. Finally, experimental results are given and analyzed in section V.

II. CREDIBILIST OCCUPANCY GRIDS

A. Grid representation

An occupancy grid is a representation of the environment which proposes to split the space into a set of cells that may be free or occupied. Usual approaches attempt to determine the probability that every cell is occupied from sensor data. The grid framework was proposed by Elfes [START_REF] Elfes | Using occupancy grids for mobile robot perception and navigation[END_REF], but it remained little used at that time because this approach is time, memory and computing consuming. During the recent years, this representation has become more and more used because it is an efficient framework for vehicle navigation since it represents in an explicit way the navigable space. Moreover, embedded computers have more powerful calculation capabilities (memory, multi-core, GPU, etc...).

Bayesian inference is the common background used to cope with errors and uncertainty. Many extensions have been published in the literature, like the Bayesian Occupancy Filter (BOF) [START_REF] Coue | Bayesian occupancy filteing for multitarget tracking : an automotive application[END_REF] which estimates simultaneously the occupation and the speed of the cells. Some works add prior information to the BOF in order to model mobile object motions [START_REF] Gindele | Bayesian occupancy grid filter for dynamic environments using prior map knowledge[END_REF]. The grid framework has currently been used to solve the complete perception problem from a mobile platform point of view in a dynamic environment [START_REF] Gate | Reliable perception of highly changing environments : Implementations for car-to-pedestrian collision avoidance systems[END_REF], since it is able to represent explicitly the free space and the locations of the obstacles.

In the present work, two grids are managed simultaneously: a ScanGrid (SG) which captures the current sensor information and a MapGrid (MG) performing the temporal integration of data in a fixed frame.

B. A credibilist framework to combine sensor data

In this work, we propose to use a Credibilist approach based on Dempster-Shafer (DS) theory because this framework offers an interesting solution to make the difference between unknown (there is no information) and doubt caused by conflicting information gathered incrementally in the fusion process. Similar ideas have been developed in [START_REF] Pagac | An evidential approach to map-building for autonomous vehicles[END_REF] to implement a DS (also called Evidential) occupancy grid intended to manage in the representation of the environment the specific uncertainty arising from the ultrasonic sensor used. Unfortunately, the DS fusion operator used in that work was very time consuming and showed a important sensitivity to aberrant measurements.

Let consider a concrete case to illustrate these concepts. There is an undiscovered cell in the map and a cell that has been observed free previously and that is now occupied by a moving object. When using a Bayesian framework, these two cells will have the same belief after the measurement update, e.g. an occupancy probability equal to 0.5. In a credibilist framework undiscovered cells receive a belief mass on a special state called "unknown" whereas in the other cases (appearance or disappearance of an object) lead to produce a belief mass on what is called "conflict".

Let see how these concepts are handled. In a way similar to probabilistic occupancy grids, every grid cell is assigned a state between two possible values F ree (F ) and Occupied (O). Let defined Ω the ignorance state in the same way that in [START_REF] Pagac | An evidential approach to map-building for autonomous vehicles[END_REF], e.g. Ω = {F, O}. In the theory of belief functions, one has to increase this set, by considering the power-set which is defined as 2 Ω = {F, O, Ω, ∅}. For each cell, a mass function is calculated and provides four beliefs on the state of the cell

[m (F ) m (O) m (Ω) m (∅)],
where m (A) represents respectively the piece of evidence that the space is free, occupied, unknown or resulting of conflict. Mass functions verify the property A⊆Ω m (A) = 1.

The Transferable Belief Model (TBM) framework proposed by Smets [START_REF] Smets | The transferable belief model[END_REF] provides also powerful tools since many fusion operators were developed in order to cope with different problems. The conjunctive combination rule, noted ∩ , is the one that is the most well known and treats the case of the fusion of two reliable sources. It is defined as follows. Let m 1 and m 2 be two given mass functions and let m 1 ∩ 2 be the result of their combination by ∩ . The result of the combination using conjunctive rule is computed in the following way:

∀A ⊆ Ω, m 1 ∩ 2 (A) = B∩C=A|B,C⊆Ω m 1 (B) • m 2 (C)
In the case of undiscovered cells by both sources S 1 and S 2 , (Eq. 1) gives the Basic Belief Assignments (BBA). The result of the combination (Eq. 2) shows that the mass remains on Ω.

∀i ∈ [1 2] m i (F ) = 0 m i (O) = 0 m i (Ω) = 1 m i (∅) = 0 (1)
m 1 ∩ 2 (F ) = 0 m 1 ∩ 2 (O) = 0 m 1 ∩ 2 (Ω) = 1 m 1 ∩ 2 (∅) = 0 (2)
In the case of two opposite observations of the sources S 1 and S 2 , the combination of the two BBA (Eq. 3) shows that the resulting mass is concentrated on the empty set illustrating the conflicting situation (Eq. 4).

m 1 (F ) = 0 m 1 (O) = 0.8 m 1 (Ω) = 0.2 m 1 (∅) = 0 m 2 (F ) = 0.8 m 2 (O) = 0 m 2 (Ω) = 0.2 m 2 (∅) = 0 (3) m 1 ∩ 2 (F ) = 0.16 m 1 ∩ 2 (O) = 0.16 m 1 ∩ 2 (Ω) = 0.04 m 1 ∩ 2 (∅) = 0.64 (4)
Contrary to a probabilistic approach, the difference between the two cases is significant. Several works have emphasized the importance of conflict analysis [START_REF] Smets | Analyzing the combination of conflicting belief functions[END_REF], [START_REF] Ramasso | State filtering and change detection using tbm conflict. application to human action recognition in athletics video[END_REF]. This particular point is a key issue in our approach.

III. FROM SENSOR DATA TO CREDIBILIST SCANGRID

In this section, we explain how to build a credibilist SG in polar coordinates using a multi-layer lidar scan with potentially several echoes in the same line of sight.

A. Sensor model

Sensors provide data from physical phenomena regardless of the contextual situation. To interpret such information, one has to model the sensor by making some simplifying assumptions. In this paper, a lidar sensor is used and we assume that the different laser beams are above the road i.e. the echoes are coming from potential obstacles. In order to be as close as possible to the sensor acquisition process, the SG is created in a polar frame.

The space is divided in several angular sectors as shown in Figure 1. Multiple measures in the same angular sector are possible because the sensor provides multi-echo measurement in several layers.

On one hand, each sector is considered independently of the others. This allows computing the SG row by row (1 row = 1 sector). On the other hand, each cell in a row depends on the other ones. This is due to the propagation of the laser beam and also to the multi-echoes strategy, since the state of cells in a row depends on their position compared to the echo position. Figure 2 illustrates the row affectation according to several lidar measurements in a sector.

B. ScanGrid computation

This section details the construction of one row of the SG, corresponding to the angular sector Θ = [θ -, θ + ]. We define a BBA for each cell computed from the sensor scan points. One cell is a box

R × Θ, R = [r -, r + ]. Let Z Θ = z i = r i θ i , θ i ∈ Θ, i ∈ [0, n] the set of n scan
points in the angular sector Θ. BBAs m S are performed in the following way.

for every cell {R, Θ} :

-∀i ∈ [1, n] if ∃ ri|ri ∈ R then mS (∅) = 0 mS (F ) = 0 mS (O) = 1 -λO mS (Ω) = λO -∀i ∈ [1, n] if ∄ ri| ri ∈ R then mS (∅) = 0 mS (F ) = 1 -λF mS (O) = 0 mS (Ω) = λF -∀i ∈ [1, n] if ∄ ri| ri ∈ R and r + < min (ri) then mS (∅) = 0 mS (F ) = 0 mS (O) = 0 mS (Ω) = 1
λ F and λ O are parameters belong to [0 , 1] and reflect the confidence in the measurement (0 if confident). This confidence is linked to the principle of measurement (false alarm or miss detection). Indeed, the angular separation between two points is larger than the width of a laser beam. So, there is a nonzero probability that an object is not detected as explained in Figure 3. It is important to notice that the non scanned angular sectors are fully affected with the vacuous mass m Ω (Eq. 5) in order to indicate that there is no information.

m Ω (F ) = 0 m Ω (O) = 0 m Ω (Ω) = 1 m Ω (∅) = 0 (5)
IV. GLOBAL FUSION AT THE MAPPING LEVEL Grid transformations, map updating and mobile objects detection are described in this section.

A. System overview

Our aim is to build a grid in a global world frame (Map-Grid) to map the navigable space and the static environment. The MG is a Cartesian map referenced in a world frame and composed of L-length squared cells. Each cell refers to an occupancy mass function defines on 2 Ω as described in section II-B. m Mi,j ,t is the mass function referring to the cell (i, j) of MG at time t.

The MapGrid is initialized using no prior information, e.g. with vacuous mass m Ω on every cell.

The updating mechanism is sequential. At a given time, the previous MapGrid, MG(t-1), is updated with the current ScanGrid SG(t) built from the sensor measurements. The result is a new MapGrid MG(t). The updating step allows classifying cells containing moving objects. The proposed scheme is composed of four phases as described by Figure 4.

Map Updating

Mobile Objects Detection

Sensor Model

B. Grid transformation

To be able to merge data from the SG, one has to transform the grid support to be spatially compatible with the MapGrid. This is achieved in two stages. First, the SG is converted from polar to Cartesian coordinates, then it is moved from the sensor frame into the reference frame of MapGrid. This requires the absolute pose of the vehicle at time t to be available. Both transformations are achieved by using bi-linear interpolation methods on each of piece of evidence of mass functions. In practice, the grid can be regarded as a multi-channel image, which allows the use of image processing tools. The result is a credibilist grid called ScanGrid+ (SG+) and mass function is noted m S + ,t .

C. Global map update with the ScanGrid

When a new ScanGrid+ SG+(t) is computed, we fuse it with the previous MapGrid MG(t-1). Every cell is considered independent from the others.

The fusion is performed using the Dempster combination rule, but it is achieved in two steps in order to keep the conflicting information.

We fuse with a conjunctive rule (Eq. 6) and we store the empty set mass m ′ (∅) for a further process. As the fusion is done element by element across the entire grid, we do not repeat the indices of the cells: m M,t refers to ∀i, j, m Mi,j ,t . Then the empty-set mass is used for normalization (Eq. 7).

m ′ M,t = m M,t-1 ∩ m S + ,t (6)
m M,t (A) = m ′ M,t (A) 1-m ′ M,t (∅) A = ∅ m M,t (∅) = 0 A = ∅ (7)
Figure 5 shows the behavior of the mass function of one cell belongs to the MG along the process composed of six situations.

In the first step, the MG cell, initialized with the vacuous mass, converges to a free state according to the SG+. Then during step 2, since no change occurs, it remains stable. Step 3 is the beginning of a change in the SG+ from F to O. Two contradictory pieces of information create a part of conflict, the other part of the mass remains on F . When SG+ remains on the same state, it forces the MG to change from F to O as shown by step 4. The following steps 5 and 6 do the inverse process.

D. Mobile object detection

If MG and SG+ contradicts themselves conflict occurs. A mobile object detection is done by analyzing the conflicting mass. This mass is computed by:

m ′ M,t (∅) = m M,t-1 (O) .m S + ,t (F )+m M,t-1 (F ) .m S + ,t (O) It is composed of two different parts.
• The term m M,t-1 (F ) .m S + ,t (O) represents the conflict produced by the fusion of an occupied cell of the SG+ with a free cell of the MG. This can occur when a moving object appears in the cell. • The term m M,t-1 (O) .m S + ,t (F ) represents the conflict produced by the fusion of a free cell of the SG+ with an occupied cell of the MG. This can occur when a moving object leaves the cell. It may be noted that the conflict is generated for a short time, because if an object remains in a new cell, it will gradually change the MG and will be considered as static.

V. EXPERIMENTAL RESULTS

The presented algorithm was tested on a data-set acquired with an equipped vehicle in urban conditions and was compared to another approach based on accumulation operators [START_REF] Moras | A lidar Perception Scheme for Intelligent Vehicle Navigation[END_REF]. The data-set is a 20-minutes long sequence, acquired in collaboration with the French Geographic Institute (IGN) in Paris. The vehicle (left Fig. 6) follows a reference track (Right Fig. 6) repeated 3 times. The results presented were conducted during the first loop. Fig. 6. Left: the test vehicle with the lidar sensor in front. Right: the track followed by the vehicle Two sensors were used: Lidar ALASCA XT and Applanix positioning system. To synchronize the data, time stamps in GPS time were used. The lidar was installed to have its lowest layer horizontally: if we assume a low pitch, the floor was not detected by lidar. The lidar frequency was 15 Hz with a maximum angular resolution of 0.25°in the vehicle front and 1°on the sides. The Applanix positioning data was post-processed in order to provide a precise pose. A camera was installed just below the lidar located on the front bumper. Images were not processed, they were just used for validation and visualization of the scene.

The SG reaches a distance of 100 m with a 180°angular aperture and a resolution of 0.5 m x 1°angle. MG covers an area of 800m x 700m with a resolution of 0.5m x 0.5m.

An offline implementation of the algorithm has been developed with Matlab and a real time version ( C++ ) is under development.

The coefficients λ F and λ o have been set at 0.5 and the detection threshold on the conflict for mobile objects detection at 0.1. Static environment mapping: Figure 7 draws the mapping results of a previously implemented method based on accumulation [START_REF] Moras | A lidar Perception Scheme for Intelligent Vehicle Navigation[END_REF]. Figure 8 draws the mapping results of the credibilist method in the same conditions. The mapping was performed in real traffic condition, i.e. with the presence of many cars and pedestrians. In order to evaluate the quality of the mapping, we use a map of road space provided by IGN (Fig. 7,8: gray cells). The black cells outside the gray area represent the building. They are mapped in the same way with both methods. The accumulation approach presents many occupied cells scattered on the free space. This phenomenas are due to moving objects. The credibilist framework is less affected because this approach offers a better management of the moving objects. Indeed, the conflict analysis allows the detection of moving objects without mapping them. So, the updating mechanism performs a filter rejecting moving objects and outliers measurements while taking advantage of redundancy. The accumulation scheme adds mobile objects into the map and clean them when they leave the place. If one object is not observed when it leaves a cell, this cell remains occupied.

Moving object detection: Figure 9 illustrates results of both method to detect mobile objects during a small sequence. This sequence shows two cars passing in front of the ego-vehicle from the right then a truck crossing from the left. As there is no clustering, the vehicle was manually circled with the same color on the camera frame and on the maps in order to identify them. Maps are displayed in global frame (MG) and give a bird view of the scene. The ego-vehicle is represented by the yellow triangle looking at the top left side.

First, let focus on the accumulation approach (middle plots on each step). On this sequence, it presents poor results. Indeed, for each vehicle, even if the front is partially detected as mobile (green cells), they are mostly miss classified (blue cells refer to static objects). We can also see that mobile objects were propagated into the map during such a time (for example: black cells on the truck). The results may be made better by tuning the classification threshold, but this will create false alarms on real static elements classified as mobile.

The credibilist approach (low plot on each step) shows correct classification on each vehicle. Mobile objects don't remain enough time on the cells to change the map. The main drawback of the method is the false alarms on static objects due to the map transformation and on free space due to outliers measurement. We can notice that with both methods, objects moving in masked area are directly mapped and not classified (cluster of occupied cells in top right for example). Finally, one can conclude that the analysis of the conflict improves greatly the performance.

VI. CONCLUSION

This paper has presented a perception fusion scheme based on both local and global grid interactions. The main contribution of this work has been to develop and test a credibilist fusion approach which provides the navigation space, the mobile obstacles and the static objects. Experimental results reported here illustrate the good performance of such a strategy for detecting mobile objects which is a crucial stage in perception. A credibilist framework seems to be an interesting way to achieve such kind of simultaneous mapping and moving objects localization using grids. Indeed, a particular attention can be given to the distinction of total ignorance with conflict arising from movements in the scene. Fig. 9. Result of the detection of moving objects, four snapshots of the sequence. Each snapshot is illustrated by a camera image of the scene and the local map displaying the moving objects. For the accumulation approach (middle) the blue refers to static and the green refers to moving. In the credibilist approach (bottom), green refers to m (∅) , the moving objects.

We have noticed a deep performance improvement compared to accumulation strategies that are quite difficult to tune and which often fail to detect large moving objects. Moreover, this strategy has shown to be very efficient for capturing pieces of information coming from all the lidar's echoes by explicitly giving no occupancy information between the echoes.

In its current implementation, precise localization is a prerequisite. One perspective is to analyze how our approach is robust to a dead-reckoning localization method, using for instance lidar odometry. Another perspective is to use another credibilist sensor model able to fuse each sensor measurement independently.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Polar space representation, showing how the information from a measurement can be interpreted. Green refers to free space, red to occupied and blue to unknown .

Fig. 3 .

 3 Fig. 3. Case of no detected object with lidar

Fig. 4 .

 4 Fig. 4. Grid based perception scheme for mapping and moving objects detection

Fig. 5 .

 5 Fig. 5. Simulation of the evolution of the mass functions. Top figure: simulated ScanGrid+ detection. Center figure: result of the conjunctive combination. Bottom figure: result of Dempster's rule.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Result of an accumulation mapping, black refers to occupied cells, gray refers to the road space provided by IGN

ACKNOWLEDGMENTS

This work is supported by the French ANR CityVip project. The system has been tested on sequences recorded by the test platform "STEREOPOLIS".