Jakob J Verbeek

Sam T Roweis

Nikos Vlassis

Non-linear CCA and PCA by Alignment of Local Models

We propose a non-linear Canonical Correlation Analysis (CCA) method which works by coordinating or aligning mixtures of linear models. In the same way that CCA extends the idea of PCA, our work extends recent methods for non-linear dimensionality reduction to the case where multiple embeddings of the same underlying low dimensional coordinates are observed, each lying on a different high dimensional manifold. We also show that a special case of our method, when applied to only a single manifold, reduces to the Laplacian Eigenmaps algorithm. As with previous alignment schemes, once the mixture models have been estimated, all of the parameters of our model can be estimated in closed form without local optima in the learning. Experimental results illustrate the viability of the approach as a non-linear extension of CCA.

Recent work on unsupervised nonlinear feature extraction has pursued several complementary directions. Various nonparametric spectral methods, such as Isomap[1], LLE[2],

. These methods, as well as ours,

Introduction

In this paper, we are interested in data that lies on or close to a low dimensional manifold embedded, possibly non-linearly, in a Euclidean space of much higher dimension. Data of this kind is often generated when our observations are very high dimensional but the number of underlying degrees of freedom is small. A typical example are images of an object under different conditions (e.g. pose and lighting). A simpler example is given in Fig. 1, where we have data in IR 3 which lies on a two dimensional manifold. We want to recover the structure of the data manifold, so that we can 'unroll' the data manifold and work with the data expressed in the underlying 'latent coordinates', i.e. coordinates on the manifold. Learning low dimensional latent representations may be desirable for different reasons, such as compression for storage and communication, visualization of high dimensional data, or as preprocessing for further data analysis or prediction tasks. deliver after training a functional mapping which can be used to convert previously unseen high dimensional observations into their low dimensional global coordinates. Like most of the above algorithms, our method performs non-linear feature extraction by minimizing a convex objective function whose critical points can be characterized as eigenvectors of some matrix. These algorithms are generally simple and efficient; one needs only to construct a matrix based on local feature analysis of the training data and then computes its largest or smallest eigenvectors using standard numerical methods. In contrast, methods like generative topographic mapping [START_REF] Bishop | GTM: the generative topographic mapping[END_REF] and self-organizing maps [START_REF] Kohonen | Self-organizing maps[END_REF] are prone to local optima in the objective function.

Our method is based on the same intuitions as in earlier work: the idea is to learn a mixture of latent variable density models on the original training data so that each mixture component acts as a local feature extractor. For example, we may use a mixture of factor analyzers or a mixture of principal component analyzers (PCA). After this mixture has been learned, the local feature extractors are 'coordinated' by finding, for each model, a suitable linear mapping (and offset) from its latent variable space into a single 'global' low-dimensional coordinate system. The local feature extractors together with the coordinating linear maps provide a global non-linear map from the data space to the latent space and back. Learning the mixture is driven by a density signal -we want to place models near the training points, while the post-coordination is driven by the idea that when two different models place significant weight on the same point, they should agree on its mapping into the global space.

Our algorithm, developed in the following section, builds upon recent work of coordination methods. As in [START_REF] Roweis | Global coordination of local linear models[END_REF], we use a cross-entropy between a unimodal approximation and the true posterior over global coordinates to encourage agreement. However we do not attempt to simultaneously learn the mixture model and coordinate since this causes severe problems with local minima. Instead, as in [START_REF] Teh | Automatic alignment of local representations[END_REF][START_REF] Brand | Charting a manifold[END_REF], we fix a specific mixture and then study the computations involved in coordinating its local representations. We extend the latter works as CCA extends PCA: rather than finding a projection of one set of points, we find projections for two sets of corresponding points {x n } and {y n } (x n corresponding to y n ) into a single latent space that project corresponding points in the two point sets as nearby as possible.

In this setting we begin by showing, in Section 3, how Laplacian Eigenmaps [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF] are a special case of the algorithms presented here when they are applied to only a single manifold. We go on, in Section 4, to extend our algorithm to a setting in which multiple different observation spaces are available, each one related to the same underlying global space but through different nonlinear embeddings. This naturally gives rise to a nonlinear version of weighted Canonical Correlation Analysis (CCA). We present results of several experiments in the same section and we conclude the paper with a general discussion in Section 5.

Non-linear PCA by aligning local feature extractors

Consider a given data set X = {x 1 , . . . , x N } and a collection of k local feature extractors, f s (x) is a vector containing the, zero or more, features produced by model s. Each feature extractor also provides an "activity signal", a s (x) representing its confidence in modeling the point. We convert these activities into posterior responsibilities using a simple softmax: p(s|x) = exp(a s (x))/ r exp(a r (x)). If the experts are actually components of a mixture, then setting the activities to the logarithm of the posteriors under the mixture will recover exactly the same posteriors above.

Next, we consider the relationship between the given representation of the data and the representation of the data in a global latent space, which we would like to find. Throughout, we will use g to denote latent 'Global' coordinates for data. For the unobserved latent coordinate g corresponding to a data point x n and conditioned on s, we assume the density:

p(g|x n , s) = N (g; κ s + A s f s (x n ), σ 2 I) = N (g; g ns , σ 2 I), (1)
where N (g; µ, Σ) is a Gaussian distribution on g with mean µ and covariance Σ. The mean, g ns , of p(g|x n , s) is the sum of the component offset κ s in the latent space and a linear transformation, implemented by A s , of f s (x n ). From now on we will use homogeneous coordinates and write:

L s = [A s κ s ] and z ns = [f s (x n ) 1]
, and thus g ns = L s z ns . Consider the posterior distribution on latent coordinates given some data:

p(g|x) = s p(s, g|x) = s p(s|x)p(g|x, s). (2)
Given a fixed set of local feature extractors and a corresponding activities, we are interested in finding linear maps L s that give rise to 'consistent' projections of the data in the latent space. By 'consistent', we mean that the p(g|x, s) are similar for components with large posterior. If the predictions are in perfect agreement for a point x n , then all the g ns are equal and the posterior p(g|x) is Gaussian, in general p(g|x) is a mixture of Gaussians. To measure the consistency, we define the following error function:

Φ({L 1 , . . . , L k }) = min {Qn,...QN } n,s q ns D(Q n (g) p(g|x n , s)), (3)
where we used q ns as a shorthand for p(s|x n ) and Q n is a Gaussian with mean g n and covariance matrix Σ n . The objective sums for each data point x n and model s the Kullback-Leibler divergence D between a single Gaussian Q n (g) and the component densities p(g|x, s), weighted by the posterior p(s|x n ). It is easy to derive that in order to minimize the objective Φ w.r.t. g n and Σ n we obtain:

g n = s q ns g ns and Σ n = σ 2 I, (4)
where I denotes the identity matrix. Skipping some additive and multiplicative constants with respect to the linear maps L s , the objective Φ then simplifies to:

Φ = n,s q ns g n -g ns 2 = 1 2 n,s,t q ns q nt g nt -g ns 2 ≥ 0. (5)
The main attraction with this setup is that our objective is a quadratic function of the linear maps L s , as in [START_REF] Teh | Automatic alignment of local representations[END_REF][START_REF] Brand | Charting a manifold[END_REF]. Using some extra notation, we obtain a clearer form of the objective as a function of the linear maps. Let:

u n = [q n1 z n1 . . . q nk z nk ], U = [u 1 . . . u N ] , L = [L 1 . . . L k ] . (6)
Note that from ( 4) and ( 6) we have: g n = (u n L) . The expected projection coordinates can thus be computed as: G = [g 1 . . . g N ] = UL. We define the block-diagonal matrix D with k blocks given by D s = n q ns z ns z ns . The objective can now be written as:

Φ = Tr{L (D -U U)L}. (7)
The objective function is invariant to translation and rotation of the global latent space and re-scaling the latent space changes the objective monotonically, c.f. (5). To make solutions unique with respect to translation, rotation and scaling, we impose two constraints:

(transl.) : ḡ = n g n /N = 0, (rot. + scale) : Σ g = n (g n -ḡ)(g n -ḡ) /N = I.
The columns of L minimizing Φ are characterized as the generalized eigenvectors:

(D -U U)v = λU Uv ⇔ Dv = (λ + 1)U Uv. ( 8
)
The value of the objective function is given by the sum of the corresponding eigenvalues λ. The smallest eigenvalue is always zero, corresponding to mapping all data into the same latent coordinate. This embedding is uninformative since it is constant, therefore we select the eigenvectors corresponding to the second up to the (d + 1) st smallest eigenvalues to obtain the best embedding in d dimensions. Note that, as mentioned in [START_REF] Teh | Automatic alignment of local representations[END_REF], this framework enables us to use feature extractors that provide different numbers of features.

In Fig. 1 we give an illustration of applying the above procedure to a simple manifold. The plots show the original data presented to the algorithm (left) and the 2-dimensional latent coordinates g n = s q ns g ns found by the algorithm (right).

Laplacian Eigenmaps as a special case

Consider the special case of the algorithm of Section 2, where no features are extracted. The only information the mixture model provides are the posterior probabilities collected in the matrix Q with [Q] ns = q ns = p(s|x n ). In that case:

g ns = κ s , U = Q, L = [κ 1 . . . κ k ] , (9)
Φ = Tr{L (D -A)L} = s,t κ s -κ t 2 n q ns q nt , (10)
where A = Q Q is an adjacency matrix with [A] st = n q ns q nt and D is the diagonal degree matrix of A with [D] ss = t A st = n q ns . Optimization under the constrains of zero mean and identity covariance leads to the generalized eigenproblem:

(D -A)v = λAv ⇔ (D -A)v = λ 1 + λ Dv (11)
The optimization problem is exactly the Laplacian Eigenmaps algorithm [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF], but applied on the mixture components instead of the data points. Since we do not use any feature extractors in this setting, it can be applied to mixture models that model data for which it is hard to design feature extractors, e.g. data that has (both numerical and) categorical features. Thus, we can use mixture densities without latent variables, e.g. mixtures of multinomials, mixtures of Hidden Markov Models, etc. Notice that in this manner the mixture model not only provides a soft grouping of the data through the posteriors, but also an adjacency matrix between the groups.

Non-linear CCA by aligning local feature extractors

Canonical Correlation Analysis (CCA) is a data analysis method that finds correspondences between two or more sets of measurements. The data are provided in tuples of corresponding measurements in the different spaces. The sets of measurements can be obtained by employing different sensors to make measurements of some phenomenon. Our main interest in this paper is to develop a nonlinear extension of CCA which works when the different measurements come from separate nonlinear manifolds that share an underlying global coordinate system. Non-linear CCA can be trained to find a shared low dimensional embedding for both manifolds, exploiting the pairwise correspondence provided by the data set. Such models can then be used for different purposes, like sensor fusion, denoising, filling in missing data, or predicting a measurement in one space given a measurement in the other space. Another important aspect of this learning setup is that the use of multiple sensors might also function as regularization helping to avoid overfitting, c.f. [START_REF] Ham | Learning high dimensional correspondences from low dimensional manifolds[END_REF].

In CCA two (zero mean) sets of points are given: X = {x 1 , . . . , x N } ⊂ IR p and Y = {y 1 , . . . , y N } ⊂ IR q . The aim is to find linear maps a and b, that map members of X and Y respectively on the real line, such that the correlation between the linearly transformed variables is maximized. This is easily shown to be equivalent to minimizing:

E = 1 2 n [ax n -by n ] 2 (12)
under the constraint that a

[ n x n x n ]a + b[ n y n y n ]b = 1.
The above is easily generalized such that the sets do not need to be zero mean and allowing a translation as well. We can also generalize by mapping to IR d instead of the real line, and then requiring the sum of the covariance matrices of the projections to be identity. CCA can also be readily extended to take into account more than two point sets, as we now show.

In the generalized CCA setting with multiple point-sets, allowing translations and linear mappings to IR d , the objective is to minimize the squared distance between all pairs of projections under the same constraint as above. We denote the projection of the n-th point in the s-th point-set as g ns and let g n = 1 k s g ns . We then minimize the error function:

Φ CCA = 1 2k 2 n,s,t g ns -g nt 2 = 1 k n,s g ns -g n 2 . ( 13
)
The objective Φ in equation ( 5) coincides with Φ CCA if q ns = 1/k. The different constraints imposed upon the optimization by CCA and our objective of the previous sections are equivalent. We can thus regard the alignment procedure as a weighted form of CCA. This suggests using the coordination technique for non-linear CCA. This is achieved quite easily, without modifying the objective function [START_REF] Bregler | Surface learning with applications to lipreading[END_REF]. We consider different point sets, each having a mixture of locally valid linear projections into the 'global' latent space that is now shared by all mixture components and point sets. We minimize the weighted sum of the squared distances between all pairs of projections, i.e. we have pairs of projections due to the same point set and also pairs that combine projections from different point sets.

We use c as an index ranging over the C different observation spaces, and write q c ns for the posterior on component s for observation n in observation space c. Similarly, we use g c ns to denote the projection due component s from space c. The average projection due to observation space c is then denoted by g c n = s q c ns g c ns . We use index r to range over all mixture components and observation spaces, so that q nr = 1 C p(s|x n ) if r corresponds to (c = 1, s) and q nr = 1 C p(s|y n ) if r corresponds to (c = 2, s), i.e. r ↔ (c, s). The overall average projection then becomes: g n = 1 C c g c n = r q nr g nr . The objective ( 5) can now be rewritten as:

Φ = n,r q nr g nr -g n 2 = 1 C c,n g n -g c n 2 + 1 C c,n,s q c ns g c n -g c ns 2 . ( 14
)
Observe how in (14) the objective sums between point set consistency of the projections (first summand) and within point set consistency of the projections (second summand). The above technique can also be used to get more stable results of the chart coordination procedure for a single manifold discussed in Section 2. Robustness for variation in the mixture fitting can be improved by using several sets of charts fitted to the same manifold. We can then align all these sets of charts by optimizing (14). This aligns the charts within each set and at the same time makes sure the different sets of aligned charts are aligned, providing important regularization, since now every point is modeled by several local models.

Note that if the charts and responsibilities are obtained using a mixture of PCA or factor analyzers, the local linear mappings to the latent space induce a Gaussian mixture in the latent space. This mixture can be used to compute responsibilities on components given latent coordinates. Also, for each linear map from the data to the latent space we can compute a pseudo inverse projecting back. By averaging the individual back projections with the responsibilities computed in latent space we obtain a projection from the latent space to the data space. In total, we can thus map from one observation space into another. This is how we generated the reconstructions in the experiments reported below. When using linear CCA for data that is non-linearly embedded, reconstructions will be poor since linear CCA can only map into a low dimensional linear subspace.

As an illustrative example of the non-linear CCA we used two point-sets in IR 2 . The first point-set was generated on an S-shaped curve the second point set was generated along an arc, see Fig. 2. To both point sets we added Gaussian noise and we learned a 10 component mixture model on both sets. In the rightmost panel of Fig. 2 the, clearly successfully, discovered latent coordinates are plotted against the coordinate on the generating curve. Below, we describe three more challenging experiments.

In the first experiment we use two data sets which we know to share the same underlying degrees of freedom. We use images of a face varying its gaze left-right and up-down. We cut these images in half to obtain our two sets of images. We trained the system on 1500 image halves of 40 × 20 pixels each. Both image halves were modeled with a mixture of 40 components. In Fig. 3 some generated right half images based on the left half are shown.

The second experiment concerns appearance based pose estimation of an object. One point set consists of a pixel representation of images of an object and the other point set contains the corresponding pose of the camera w.r.t. the object. For the pose parameters we used the identity to 'extract' features (i.e. we just used one component for this space). The training data was collected1 by moving a camera over the half-sphere centered at the object. A mixture of 40 PCA's was trained on the image data and aligned with the pose parameters in a 2-dimensional latent space. The right panel of Fig. 3 shows reconstructions of the images conditioned on various pose inputs (left image of each pair is reconstruction based on pose of right image). Going the other way, when we input an image and estimate the pose, the absolute errors in the longitude (0 • -360 • ) were under 10 • in over 80% of the cases and for latitude (0 • -90 • ) this was under 5 • in over 90% of the cases. In the third experiment we use the same images as in the second experiment, but replace the direct (low dimensional) supervision signal of the pose parameters with (high dimensional) correspondences in the form of images of another object in corresponding poses. We trained a mixture of 40 PCA's on both image sets (2000 images of 64×64 pixels in each set) and aligned these in a 3-dimensional latent space. Comparing the pose of an object to the pose of the nearest (in latent space) image from the other object the std. dev. of error in latitude is 2.0 • . For longitude we found 4 errors of about 180 • in our 500 test cases, the rest of the errors had std. dev. 3.9 • . Given a view of one object we can reconstruct the corresponding view of the second object, Fig. 4 shows some of the obtained reconstruction results. All presented reconstructions were made for data not included in training.

Discussion

In this paper, we have extended alignment methods for single manifold nonlinear dimensionality reduction to perform non-linear CCA using measurements from multiple manifolds. We have also shown the close relationship with Laplacian Eigenmaps [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF] in the degenerate case of a single manifold and feature extractors of zero dimensionality.

In [START_REF] Teh | Automatic alignment of local representations[END_REF] a related method to coordinate local charts is proposed, which is based on the LLE cost function as opposed to our cross-entropy term; this means that we need more than just a set of local feature extractors and their posteriors: we also need to be able to compute reconstruction weights, collected in a N × N weight matrix. The weights indicate how we can reconstruct each data point from its nearest neighbors. Computing these weights requires access to the original data directly, not just through the "interface" of the mixture model. Defining sensible weights and the 'right' number of neighbors might not be straightforward, especially for data in non-Euclidean spaces. Furthermore, computing the weights costs in principle O(N 2 ) because we need to find nearest neighbors, whereas the presented work has running time linear in the number of data points.

In [START_REF] Ham | Learning high dimensional correspondences from low dimensional manifolds[END_REF] it is considered how to find low dimensional representations for multiple point sets simultaneously, given few correspondences between the point sets. The generalization of LLE presented there for this problem is closely related to our non-linear CCA model. The work presented here can also be extended to the case where we know only for few points in one set to which points they correspond in the other set. The use of multiple sets of charts for one data set is similar in spirit as the self-correspondence technique of [START_REF] Ham | Learning high dimensional correspondences from low dimensional manifolds[END_REF] where the data is split into several overlapping sets used to stabilize the generalized LLE. Finally, it would be interesting to compare our approach with treating the data in the joint (x, y) space and employing techniques for a single point set [START_REF] Brand | Charting a manifold[END_REF][START_REF] Teh | Automatic alignment of local representations[END_REF][START_REF] Roweis | Global coordination of local linear models[END_REF]. In this case, points for which we do not have the correspondence can be treated as data with missing values.

Figure 1 :

 1 Figure 1: Data in IR 3 with local charts indicated by the axes (left). Data representation in IR 2 generated by optimizing our objective function. Expected latent coordinates g n are plotted (right).

Figure 2 :

 2 Figure 2: Data and charts, indicated by bars (left-middle). Latent coordinates (vert.) and coordinate on generating curve (hor.) (right).

Figure 3 :

 3 Figure 3: Right half of the images was generated given the left half using the trained model (left). Image reconstructions given pose parameters (right).

  Figure 4: I1 : image in first set (a), I2 : corresponding image in second set (b), closest image in second set (in latent space) to I1 (c), reconstruction of I2 given I1 (d).

Thanks to G. Peters for sharing the images used in[START_REF] Peters | How to measure the pose robustness of object views[END_REF] and recorded at the Institute for Neural Computation, Ruhr-University Bochum, Germany.

Acknowledgments

JJV and NV are supported by the Technology Foundation STW (AIF4997) applied science division of NWO and the technology program of the Dutch Ministry of Economic Affairs. STR is supported in part by the Learning Project of IRIS Canada and by the NSERC.


Jan R J Nunnink
email: jnunnink@science.uva.nl

Jakob J Verbeek
email: jverbeek@science.uva.nl

Nikos Vlassis
email: vlassis@science.uva.nl

Accelerated Greedy Mixture Learning

Mixture probability densities are popular models that are used in several data mining and machine learning applications, e.g., clustering. A standard algorithm for learning such models from data is the Expectation-Maximization (EM) algorithm. However, EM can be slow with large datasets, and therefore approximation techniques are needed. In this paper we propose a variational approximation to the greedy EM algorithm which offers speedups that are at least linear in the number of data points. Moreover, by strictly increasing a lower bound on the data log-likelihood in every learning step, our algorithm guarantees convergence. We demonstrate the proposed algorithm on a synthetic experiment where satisfactory results are obtained.

Introduction

An important task in data analysis is to search and identify clusters of data items. These clusters can be used for a number of purposes, such as classification of new data or for predicting missing data. Data clustering methods can also be used as learning algorithms when applied to partially labelled data, where each found cluster corresponds to a class.

The task of clustering algorithms is to fit a model to the data. A useful and popular class of models are Mixture models which are convex combinations of basic model components [START_REF] Mclachlan | Finite Mixture Models[END_REF]. Typically these components are Gaussian density functions, because often data is locally Gaussian distributed, in which case the data set can be assumed to have been generated by a hypothetical Gaussian mixture. The task is then to find the parameters of that generating mixture.

Because of their probabilistic nature, Gaussian mixtures are in principle preferred over models that partition a data set in discrete parts. In most applications where a new data item needs to be classified, it is more desirable to calculate the probability that this item belongs to certain clusters than to assign it to strictly one cluster.

The most popular algorithm for training a Gaussian mixture is the Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. Among its advantages are its easy implementation, no need to set extra user-defined parameters, and guaranteed monotone increase in model quality. There are also some downsides, however, the most important of which are its high dependency on initialization, and its computational complexity which is linear with respect to the size of the data set.

Over the past years several improvements of the EM algorithm were proposed to counteract those problems. One of them is a generalization of the EM algorithm, analogous to the (negative) variational free energy in statistical physics, that justifies several useful variants [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]. Also, some improvements aimed at lowering the computational cost by working with groups of data items instead of the items themselves [START_REF] Moore | Very fast EM-based mixture model clustering using multiresolution kd-trees[END_REF] (In [START_REF] Alsabti | An efficient k-means clustering algorithm[END_REF][START_REF] Kanungo | An efficient k-means clustering algorithm: Analysis and implementation[END_REF] this method is applied to the k-means algorithm). Finally, a greedy version of EM was proposed, which deals with the initialization problems of EM, thus resulting in higher quality models (Verbeek et al., 2003;Vlassis and Likas, 2002). This greedy method, however, is computationally more expensive.

In this paper we propose a variational approximation to the greedy EM algorithm for Gaussian mixtures proposed in (Verbeek et al., 2003) that offers a speedup without compromising stability. As in [START_REF] Moore | Very fast EM-based mixture model clustering using multiresolution kd-trees[END_REF], we first partition the data and cache some statistics in each partition cell. A variational bound [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] allows us to compute for each cell optimal responsibilities over mixing components (E-step), which are further used to update the mixture parameters (M-step). Both steps have cost that is independent of the size of the data set, and is linear in the number of partition cells. We demonstrate the proposed algorithm on a synthetic experiment where very satisfactory results are obtained. A more detailed treatment can be found in [START_REF] Nunnink | Large scale Gaussian mixture modelling using a greedy Expectation-Maximisation algorithm[END_REF].

Gaussian mixtures and the EM algorithm

A k-component Gaussian mixture for a random vector x in IR d is defined as the convex combination

p(x) = k s=1 p(x|s) p(s) (1) of d-dimensional Gaussian densities p(x|s) = (2π) -d/2 |C s | -1/2 exp - 1 2 (x -m s ) C -1 s (x -m s ) (2)
each parameterized by its mean m s and covariance matrix C s . The components of the mixture are indexed by the random variable s that takes values from 1 to k, and p(s) defines a discrete 'prior' distribution over the components. Given a set {x 1 , . . . , x n } of points, assumed to be independent and identically distributed, the learning task is to estimate the parameter vector θ = {p(s), m s , C s } k s=1 of the k components that maximizes the log-likelihood function

L(θ) = n i=1 log p(x i ; θ) = n i=1 log k s=1 p(x i |s)p(s).
(3) Throughout we assume that the likelihood function is bounded from above (e.g., by placing lower bounds on the eigenvalues of the components covariance matrices) in which case the maximum likelihood estimate is known to exist [START_REF] Lindsay | The geometry of mixture likelihoods: a general theory[END_REF].

The variational EM algorithm

Maximization of the data log-likelihood L(θ) can be efficiently carried out by the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF].

In this work we consider a variational generalization of EM [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] which allows for useful extensions.

The variational EM algorithm performs iterative maximization of a lower bound of the data log-likelihood. In our case, this bound F(θ, Q) is a function of the current mixture parameters θ and a factorized distribution Q = n i=1 q i (s), where each q i (s) corresponds to a data point x i and defines an arbitrary discrete distribution over s. For a particular realization of s we will refer to q i (s) as the 'responsibility' of component s for the point x i .

This lower bound, analogous to the (negative) variational free energy in statistical physics, can be expressed by the following two (equivalent) decompositions

F(θ, Q) = n i=1 [log p(x i ; θ) -D(q i (s) p(s|x i ; θ))] (4) = n i=1 k s=1 q i (s)[log p(x i , s; θ) -log q i (s)] (5)
where D(• •) denotes the Kullback-Leibler divergence between two distributions, and p(s|x i ) is the posterior distribution over components of a data point x i computed from (1) by applying Bayes' rule. The dependence of p on θ is throughout assumed, although not always written explicitly.

Since the Kullback-Leibler divergence between two distributions is non-negative, the decomposition (4) defines indeed a lower bound on the log-likelihood. Moreover, the closer the responsibilities q i (s) are to the posteriors p(s|x i ), the tighter the bound. In particular, maxima of F are also maxima of L [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]. In the original derivation of EM [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], each E step of the algorithm sets q i (s) = p(s|x i ) in which case, and for the current value θ t of the parameter vector, holds F(θ t , Q) = L(θ t ). However, other (suboptimal) assignments to the individual q i (s) are also allowed provided that F increases in each step.

For particular values of the responsibilities q i (s), we can solve for the unknown parameters of the mixture by using the second decomposition (5) of F. It is easy to see that maximizing F for the unknown parameters of a component s yields the following solutions:

p(s) = 1 n n i=1 q i (s), (6)
m s = 1 np(s) n i=1 q i (s)x i , (7)
C s = 1 np(s) n i=1 q i (s)x i x i -m s m s , (8)
from which we directly see the linear complexity of EM in the number n of datapoints.

Locally shared responsibilities

As mentioned above, in each step of the variational EM we are allowed to assign any responsibilities q i (s) to the data as long as this increases F. The idea behind our algorithm is to assign equal responsibilities to chunks of data points that are nearby in the input space.

Consider a partitioning P of the data space into a collection of non-overlapping cells {A 1 , . . . , A m }, such that each point in the data set belongs to a single cell. To all points in a cell A ∈ P we assign the same responsibility distribution q A (s) which we can compute in an optimal way as we show next.

Note from (5) that the objective function F can be written as a sum of local parts F = A∈P F A , one per cell. If we impose q i (s) = q A (s) for all data points x i ∈ A, then the part of F corresponding to a cell A reads

F A = x i ∈A k s=1 q A (s)[log p(x i |s) + log p(s)- log q A (s)] = |A| k s=1 q A (s)[log p(s) -log q A (s)+ 1 |A| x i ∈A log p(x i |s)]. (9)
If we set the derivatives of F A w.r.t. q A (s) to zero we find the optimal distribution q A (s) that (globally) maximizes F A :

q A (s) ∝ p(s) exp log p(x|s) A (10)
where • A denotes average over all points in cell A. Such an optimal distribution can be separately computed for each cell A ∈ A, and only requires computing the average joint loglikelihood of the points in A.

Speedup using cached statistics

We now show that it is possible to efficiently compute (i) the optimal q A (s) for each cell A in the E-step and (ii) the new values of the unknown mixture parameters in the M-step, if some statistics of the points in each cell A are cached in advance. The averaging operation in (10) can be written1 :

log p(x|s) A = 1 |A| x i ∈A log p(x i |s) = - 1 2 [log |C s | + 1 |A| x i ∈A (x i -m s ) C -1 s (x i -m s )] = - 1 2 [log |C s | + m s C -1 s m s + x C -1 s x A - 2m s C -1 s x A ] = - 1 2 [log |C s | + m s C -1 s m s + Tr{C -1 s xx A }- 2m s C -1 s x A ]. ( 11
)
¿From this we see that the mean x A and covariance xx A of the points in A are sufficient statistics for computing the optimal responsibilities q A (s).

The same statistics can be used for updating the mixture parameters θ. If we set the derivatives of (9) w.r.t. θ to zero we find the update equations:

p(s) = 1 n A∈P |A|q A (s), (12)
m s = 1 np(s) A∈P |A|q A (s) x A , (13)
C s = 1 np(s) A∈P |A|q A (s) xx A -m s m s , (14)
in direct analogy to the update equations ( 6)-( 8), with the advantage that the linear complexity in the number of data points has been replaced by linear complexity in the number of cells of the particular partitioning. Note that, whatever partitioning we choose, the Chunky EM algorithm presented above (which interacts with the data only through the cached statistics of chunks of data) is always a convergent algorithm that strictly increases in each step a lower bound on data log-likelihood. In the limit, if we partition all data points into separate cells, the algorithm will converge to a local maximum of the data log-likelihood.

Greedy mixture learning

A recent approach to mixture learning involves building a mixture in a 'greedy' manner (Verbeek et al., 2003;Vlassis and Likas, 2002;[START_REF] Li | Mixture density estimation[END_REF]. The idea is to start with a single component (which is trivial to find), and then alternate between adding a new component to the mixture and updating the complete mixture.

In particular, given a k-component Gaussian mixture p k (x) that has converged, the greedy method seeks a new component φ(x) with mean m φ and covariance C φ , and a mixing weight a ∈ (0, 1) that maximize the log-likelihood L k+1 of the two-component mixture p k+1 (x) = (1 -a)p k (x) + aφ(x; m φ , C φ ), (15) where p k (x) is kept fixed. The advantages of greedy mixture learning are: (1) initializing the mixture is trivial, (2) local maxima of L are easier to escape, and (3) model selection becomes more manageable.

In (Verbeek et al., 2003), the search for a good component to add to p k (x) involves first splitting the data according to their 'nearest' (with highest posterior) component, then randomly generating a number of candidate components from the points in each subset, and finally maximizing L k+1 using only the data from the corresponding subset.

The same principle can also be applied in the case of pre-partitioned data sets. In particular, in component allocation we divide all cells A ∈ P into k disjoint subsets P s according to their 'nearest' (with highest responsibility) com-ponent:

P s = {A ∈ P : s = arg max s q A (s )}.
(16)

Then we randomly split P s in two sets of adjacent cells and for each such set S we generate a component φ(x; m φ , C φ ) from the data contained in S:

m φ = 1 n S A∈S n A x A , (17)
C φ = 1 n S A∈S n A xx A -mm , (18)
where n S is the total number of points in S.

Note that both m φ and C φ can be calculated without requiring the data points themselves. Subsequently we update (a, m φ , C φ ) in ( 15) by maximizing a lower bound of L k+1 using only the cells in P s . (Cells outside P s will not contribute significantly to the bound.) Let q A be the responsibility of the new component φ(x; m φ , C φ ) for any cell A ∈ P s , and 1 -q A the responsibility of the old mixture p k . The free energy (9) for cell A under the two-component mixture ( 15) then reads

F k+1 A = n A q A log a q A + log φ(x) A + n A (1 -q A ) log 1 -a 1 -q A + log p k (x) A . ( 19
)
Since F k A is a lower bound on x∈A log p k (x) which we have already computed from (9), we can replace the latter in ( 19) to get the bound

F k+1 A ≥ n A q A log a q A + log φ(x) A + n A (1 -q A ) log 1 -a 1 -q A + F k A n A . (20)
In the E-step we compute the optimal q A for each cell A ∈ P s by setting the derivative of (20) w.r.t. q A to zero. This gives:

q A = a exp log φ(x) A (1 -a) exp(F k A /n A ) + a exp log φ(x) A , (21
) where log φ(x) A can be computed fast using (11). Similarly, in the M-step we maximize F k+1 P = A∈P F k+1 A using the q A found in (21).

As in (Verbeek et al., 2003), we set the responsibility of the new component for all cells outside P s to zero, in which case it is not difficult to see that we get the following update equations:

a = A∈Ps n A q A n , (22)
m φ = A∈Ps n A q A x A na , (23)
C φ = A∈Ps n A q A xx A na -m φ m φ . ( 24
)
Note that the sums run over cells in P s ⊂ P.

We refer to [START_REF] Nunnink | Large scale Gaussian mixture modelling using a greedy Expectation-Maximisation algorithm[END_REF] for more details.

Choosing a partition

The analysis presented in the previous sections applies to any partition, as long as sufficient statistics of the data have been stored in the corresponding cells. As we showed above, for any partition we obtain a convergent algorithm that strictly increases in each step a lower bound on the data log-likelihood. Moreover, by refining a given partition F cannot decrease, and in the limit holds F = L making the approximation bound tight. Clearly, various trade-offs can be made between the computational cost and the approximation quality.

A convenient structure for storing statistics in a way that permits the use of different partitions in the course of the algorithm is a kdtree [START_REF] Bentley | Multidimensional binary search trees used for associative searching[END_REF][START_REF] Moore | Very fast EM-based mixture model clustering using multiresolution kd-trees[END_REF]. This is a binary tree in which the root contains all data points, and each node is recursively split by a hyperplane that cuts through the data points contained in the node. Typically, axis-aligned hyperplanes are used for splitting nodes. In our experiments we used hyperplanes that cut along the bisector of the first principal component of the points in the node, leading to irregularly shaped cells (Sproull, 1991). As in [START_REF] Moore | Very fast EM-based mixture model clustering using multiresolution kd-trees[END_REF], we store in each node of the kd-tree the sufficient statistics of all data points under this node. Building the kd-tree and storing statistics in its nodes has cost O(n log n), but this needs to be done only once at the beginning of the algorithm.

The outer nodes of a given expansion of the kd-tree form a partition P of the data set. Further expanding the tree means refining a current partition. In our implementations as heuristic to guide the tree expansion we employ a bestfirst search strategy in which we expand the node that leads to maximal increase in F. Note that computing the change in F involves only a node and its children so it can be done fast. We also need a criterion when to stop expanding the tree, and one could use among others bounds on the variation of the data posteriors inside a node like in [START_REF] Moore | Very fast EM-based mixture model clustering using multiresolution kd-trees[END_REF], a bound on the size of the partition (number of outer nodes at any step), or sampled approximations of the difference between log-likelihood and F. Another possibility, which we adopted in our experiments, is to control the tree expansion based on the performance of the algorithm, that is, we refine a partition only if this (significantly) improves the value of F.

Related work

The idea of using a kd-tree for accelerating the EM algorithm in large-scale mixture modeling was first proposed in [START_REF] Moore | Very fast EM-based mixture model clustering using multiresolution kd-trees[END_REF]. In that work, in each EM step every node in the kdtree is assigned responsibility distribution equal to the Bayes posterior of the centroid of the data points stored in the node, i.e., q A = p(s| x A ). (Compare this with (10).) If there is little variation in the posteriors within a node, which is in turn achieved by having relatively fine partitions, the approximation q A = p(s| x A ) will only slightly affect the update in the M-step and therefore this will probably increase the data log-likelihood. However this is not guaranteed. Also, in [START_REF] Moore | Very fast EM-based mixture model clustering using multiresolution kd-trees[END_REF] a different tree expansion is computed in each EM step, while as stopping criterion for tree expansion bounds are used on the variation of the posterior probabilities of the data inside a node of the kd-tree (a nontrivial operation that in principle would require solving a quadratic programming problem).

The main advantage of our method compared to [START_REF] Moore | Very fast EM-based mixture model clustering using multiresolution kd-trees[END_REF] is that our algorithm strictly increases in each step a lower bound of the data log-likelihood by computing the optimal responsibility distribution for each node. Moreover, this optimal distribution is independent of the size, shape, or other properties of the node, allowing us to use even rough partitions. As mentioned above and as demonstrated in the experiments below, by gradually refining the partition while running the algorithm we can get close to the optima of the log-likelihood in relatively few steps.

Demonstration

In a synthetic experiment we compared the greedy EM algorithm described in (Verbeek et al., 2003) with the proposed accelerated greedy algorithm, for learning a k-component Gaussian mixture. We used artificially generated data sampled from a randomly initialized kcomponent Gaussian mixture in d dimensions with a component separation of c [START_REF] Dasgupta | Learning mixtures of Gaussians[END_REF]. For each data set we built a kd-tree and stored in its nodes the data statistics as explained above.

We started the accelerated greedy EM algorithm with a partition of size four and expanded the tree one node at a time, best first, as in the first experiment. We used a default data set of 10,000 points and a test set of 500 points drawn from a 5-component 2-separated mixture in two dimensions. In Fig. 1 we show the results averaged over 20 trials. The accelerated greedy algorithm is always faster, with a speedup that is linear in the size of the data set. Moreover, this speedup comes almost 'for free': the loglikelihoods of both algorithms are practically equal to that of the generating mixture.

Conclusions

We have presented a variant of the Expectation-Maximization algorithm specifically designed for dealing with large-scale data sets, without paying much in terms of quality. The experiment described above clearly shows that the proposed algorithm is faster than previous versions and that this speedup is linear with respect to the size of the data set. The experiment also shows that there is barely any loss in quality even when working with large or complex data sets.

The theory shows that the effect of the EM algorithm, namely the guaranteed monotone increase in model quality with respect to the free energy, still holds independent of the size of the partition. Therefore one can use the size of the partition as an easy to implement trade-off between speed and quality. This is useful when there are, for example, only limited resources available.

Figure 1 :

 1 Figure1: Greedy mixture learning with the algorithm in(Verbeek et al., 2003) vs. the accelerated greedy EM. The graphs show the speedup factor and the bar charts show the negative log-likelihood: generating mixture (black), greedy EM (light), accelerated greedy EM (dark).

We ignore the additive constantd

log(2π) which translates into a multiplicative constant in (10).

Acknowledgements

This research is supported by the Technology Foundation STW (project AIF 4997) applied science division of NWO and the technology program of the Dutch Ministry of Economic Affairs.


Nikos Vlassis
email: vlassis@science.uva.nl

Jakob J Verbeek

intelligent autonomous systems

Keywords: Gaussian mixture, EM algorithm, bound optimization, noisy data

We address the problem of learning a Gaussian mixture from a set of noisy data points. Each input point has an associated covariance matrix that can be interpreted as the uncertainty by which this point was observed. We derive an EM algorithm that learns a Gaussian mixture that minimizes the Kullback-Leibler divergence to a variable kernel density estimator on the input data. The proposed algorithm performs iterative optimization of a strict bound on the Kullback-Leibler divergence, and is provably convergent.

The setup

We are given a set {x 1 , . . . , x n } of measured data points in IR d together with their uncertainties {C 1 , . . . , C n }, where C i is a d × d covariance matrix associated with measurement x i . The task is to learn (fit the parameters of) a k-component Gaussian mixture p(x) from the data, where

p(x) = k s=1 p(x|s)p(s), (1)
and where p(x|s) is the d-variate Gaussian density

p(x|s) = |S s | -1/2 (2π) d/2 exp - 1 2 (x -m s ) S -1 s (x -m s ) , (2)
parameterized by its mean m s and covariance matrix S s . The components of the mixture are indexed by the random variable s that takes values from 1 to k, and p(s) defines a discrete prior distribution over the components. The set of parameters we want to estimate will be denoted θ = {p(s), m s , S s } k s=1 . Given that each data point has a covariance associated with it, we can define a variable kernel density estimator [START_REF] Terrell | Variable kernel density estimation[END_REF] in the form

f (x) = 1 n n j=1 f (x|j), (3)
where f (x|j) is a d-variate Gaussian like in ( 2) with known mean x j and covariance C j , and where the index j runs over the complete data set.

The learning problem

The learning problem can be expressed as an optimization problem where we are looking for the parameter vector that minimizes the Kullback-Leibler divergence between the kernel estimate and the unknown mixture, that is

θ * = arg min θ D[f (x) p(x; θ)], (4)
where the Kullback-Leibler divergence reads

D[f (x) p(x)] = x dxf (x) log f (x) p(x) = const - x dxf (x) log p(x) (5)
and where we used the fact that f is fixed. Dropping constants, the objective function to maximize reads

L = n j=1 x dxf (x|j) log p(x) (6)
which can be regarded as an average log-likelihood function over the kernels f (x|j).

An EM approach

In order to maximize L we consider an EM approach in which we iteratively maximize a lower bound of L [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]. This bound F is a function of the current mixture parameters θ and a set of distributions (or 'responsibilities') q j (s), for j = 1, . . . , n, and it is computed by subtracting from each log-likelihood term a positive quantity (a KL-divergence) as follows:

F = n j=1 x dxf (x|j) log p(x) -D[q j (s) p(s|x)] (7)
which can be rewritten as

F = n j=1 k s=1 q j (s)
x dxf (x|j) log p(x|s) + log p(s) -log q j (s) .

The integral can be analytically computed (dropping constants) as

x dxf (x|j) log p(x|s) = - 1 2 log |S s | + (x j -m s ) S -1 s (x j -m s ) + Tr[S -1 s C j ] (9)
where Tr denotes the trace of a matrix.

The update equations

We can now directly maximize F w.r.t. the responsibilities q j (s) and the unknown mixture parameters θ. The resulting update equations are very similar to those of the EM algorithm for noise-free data. For the responsibilities we have

q j (s) ∝ p(s|j) exp - 1 2 Tr[S -1 s C j ] (10)
where p(s|j) is the Bayes posterior as in the standard EM:

p(s|j) = p(x j |s)p(s) p(x j ) . ( 11
)
For the mixture parameters we get

p(s) = n j=1 q j (s) n , m s = n j=1 q j (s)x j np(s) , S s = n j=1 q j (s)(x j x j + C j ) np(s) -m s m s . ( 12
)
The terms x j x j + C j can be computed in advance. Note that mixture learning with noisy data automatically achieves a sort of regularization: the matrix S s will always be non-singular as long as the individual matrices C j are non-singular.


R Rouveure

P Faure

M O Monod

Radar-based SLAM without odometric sensor

For environmental management, remote sensing techniques are of major interest as they can provide global and local coverage of natural environments. Considering the visibility conditions that are often poor in outdoor environments (day/night cycles, clouds, shadows, fog, etc.), a microwave system offers an alternative solution in order to overcome the difficulties of use of optical sensors. In order to deal with the need of high spatial resolution, close temporal acquisition, and low-cost data availability, we developed a small-sized microwave radar sensor called K2Pi. This radar is based on the principle of Frequency-Modulated Continuous Wave (FM-CW) technology. An application in simultaneous localization and mapping using scan matching of radar images (R-SLAM algorithm) has been developed in order to build radar map of the environment. Avoiding the use of odometric sensors, this SLAM algorithm allows placing the radar on every kind of vehicles. We present in this paper its potentialities through an application of natural environment mapping.

Introduction

Anthropogenic or natural processes have important effects on natural environments evolution. The monitoring of these evolutions is crucial for global change research (large scale) [START_REF] Nemani | Satellite monitoring of global land cover changes and their impact on climate[END_REF], but is also important at a local scale because of direct and imminent influences on vulnerable and small ecosystems or human life [START_REF] Lambin | Land-Use and Land-Cover Change -Local Processes and Global Impacts[END_REF]). In the domain of sustainable development, the environment risk diagnosis or the monitoring of particular environmental parameters imply to develop environmental models in order to estimate and predict the consequences of specific human behaviors or natural phenomena. The development of numerical models implies the availability of accurate and redundant data. But robust solutions for perception in extended outdoor environments are not so frequent, particularly with regard to the problem of resistance to the environmental conditions [START_REF] Hague | Ground based sensing systems for autonomous agricultural vehicles[END_REF][START_REF] Hellström | Autonomous Navigation for Forest Machines[END_REF]). Laser and video based systems are common optical based sensors, particularly well adapted in indoor environments. However, these systems are affected by visibility conditions that are often poor in outdoor environments: day/night cycles change illumination, weather phenomena as fog, rain or dust may deteriorate the quality of information provided. In the microwave range, data can be acquired independently of atmospheric conditions and time of the day. So radar can provide an alternative solution to overcome the shortcomings of optical sensors. Moreover, backscattering signal has an unquestionable potential for environmental applications, as it can be shown in satellite remote sensing [START_REF] Hallikainen | A Helicopter-Borne Eight-Channel Ranging Scatterometer for Remote Sensing: Part I: System Description[END_REF]. In response to new requirements related to the monitoring of environment at lower scales in comparison to airborne/spaceborne radar systems, we developed a small-sized radar sensor easy to embark on board of different vehicles and we show in this paper its potentialities through a preliminary application in natural environment mapping. Because of its intrinsic characteristics, microwave radar is particularly well suited for outdoor applications and large scale mapping, and particularly in outdoor mobile robotics. Microwave radar has been successfully used in numerous domains such as terrain imaging [START_REF] Scheding | Terrain imaging and perception using millimetre wave radar[END_REF][START_REF] Brooker | Using Millimetre Wave Radar in Underground Mines[END_REF][START_REF] Brooker | A W-band airborne interrupted FMCW imaging radar[END_REF], target tracking [START_REF] Brooker | Millimetre Wave Radar for Tracking and Imaging applications[END_REF], three-dimensional map building [START_REF] Foessel | Three dimensional map building with mmw radar[END_REF], vehicle localization and vehicle navigation [START_REF] Clark | The design of a high performance mmw radar system for autonomous land vehicle navigation[END_REF]. Moreover in recent years different simultaneous localization and mapping algorithms based on microwave radar sensor have been developed [START_REF] Dissanayake | A solution to the simultaneous localisation and map building (SLAM) problem[END_REF][START_REF] Jose | Millimetre Wave RADAR Spectra Simulation and Interpretation for Outdoor SLAM[END_REF][START_REF] Foessel | Radar Sensor Model for Three-Dimensional Map Building[END_REF]. The principle and the characteristics of the K2Pi radar imager is presented in Section 2. Section 3 develops and describes the R-SLAM algorithm which is based on a scan matching approach and the description of the environment with an occupancy grid. Two results obtained in natural environments are presented in Section 4.

Presentation of the K2Pi radar imager

Radar principle

For short-range distance measurements i.e. few hundreds of meters, the classical pulse radar technology is not the best adapted solution. For example, a precision of 1 meter in distance measurement imposes a precision of 3.3 ns in time measurement, therefore fast and expensive electronics. We developed the K2Pi microwave radar (meaning 360° scanning in K band) using the frequency modulation continuous wave (FMCW) principle. This principle is known and used for several decades [START_REF] Collette | Un altimètre radioélectrique à modulation de fréquence[END_REF][START_REF] Monod | Frequency modulated radar: a new sensor for natural environment and mobile robotics[END_REF]. When the range is some hundreds of meters, frequency modulation presents the following advantages: -low transmission power, as it is the mean power which determines the range, which is safer for the user, -transposition of the temporal variables into the frequency domain: a very short delay time ∆t is transformed into a broad variation of frequency ∆f, easier to measure.

The principle of FMCW radar consists to mix the emitted signal following a linear modulation waveform, with the received signal delayed after reflecting on the targets present within the view of the antenna [START_REF] Skolnik | Introduction to radar systems[END_REF]. Let us consider a sawtooth modulation at frequency F m over a range ∆F and a distribution of i stationary targets at respective distances R i in the antenna beam. When the echo signal is heterodyned with a portion of the transmitted signal, a beat signal S b will be produced [START_REF] Monod | Frequency modulated radar: a new sensor for natural environment and mobile robotics[END_REF][START_REF] Rouveure | Multiple targets detection with a FMCW radar dedicated to mobile robotics[END_REF]):

∑         Φ +       ∆ = i i i m ri e b t c R F F V V k S 2 2 cos π . ( 1
)
V e is the amplitude of emitted signal, V ri and Φ i respectively the amplitude of received signal and a phase term depending on target i, c the light velocity and k a mixer coefficient. When more than one target is present within the view of the radar, the mixer output will contain a sum ΣF bi of frequency components, each of them corresponding to a particular target i:

c R F F F i m bi ∆ = 2 . ( 2
)
The range R i to each target may be determined by measuring the individual frequency components by the mean of a Fourier transform. At each angle of sight a power spectrum is computed. Contrary to laser measurement which gives the nearest echo received in the beam, the aperture of the radar beam allows multi-targets detection in the same angle-ofsight. An example of radar power spectrum is presented in Fig. 1. We will notice that moving targets which introduce a Doppler frequency shift are not considered in this paper.

Radar image construction

A so-called "panoramic radar image" covering 360° i n the horizontal plane is carried out using a PPI (Plan Position Indicator) representation, which displays range of targets versus angle. Considering the combined effect of antenna rotation and vehicle motion (translation and rotation), one would observe intra-scan distortions on the final 360° image. To achieve a vehicle motion correction, two proprioceptive sensors are used to estimate the displacement of the vehicle and to take into account these distortions: an odometer to measure the longitudinal displacement of the vehicle and a gyrometer to measure the rotation of the vehicle.

Once the vehicle motion has been corrected, an anti-speckle filter is added. The speckle effect is a random amplitude modulation which exists in any radar image construction. Speckle noise reduction can be achieved either by spatial filtering or multi-look filtering [START_REF] Lee | Speckle filtering of synthetic aperture radar images: a review[END_REF][START_REF] Guoqing | Study on speckle reduction in multi-look polarimetric SAR image[END_REF]). We chose the latter approach: n independent images of the scene are incoherently combined (by space-domain averaging) to produce a single final image.

Lack of odometric sensor

Depending on the vehicle or the environment, odometric sensors can be unavailable. For example if we consider a boat on a river, it is placed on a "moving environment" and no odometric sensor is available to estimate the displacement of the radar (odometric information are used during the Vehicle Motion Correction and the Anti-Speckle Filter steps). Data from a kinematic GPS system (with a centimetric precision) could be employed, but we try to avoid the use of such a system for different reasons:

-we want to build a radar imager independent of any localization system, -the cost of a kinematic GPS system is very high, -if no differential GPS correction signal is available, it is necessary to install its own GPS ground station and the cost of the GPS system is doubled, -considering the GPS satellite constellation, the terrain topography, the presence of high vegetation or various obstacles, etc., GPS signals can be intermittently available or completely unavailable. So a direct estimation process of the vehicle displacement was developed: it allows estimating of the movement of the vehicle within initial radar images (before vehicle motion and anti-speckle corrections).

The K2Pi radar is monostatic, i.e. a common antenna is used for both transmitting and receiving. The rotating antenna achieves a complete 360° scan around the vehicle in one second, and a radar acquisition is realized each degree. The sensor includes microwave components, electronic devices for emission and reception, data acquisition and signal processing units. A general view of the K2PI radar is shown in Fig. 2. Data acquisition and signal processing units are based on a Pentium M / 2 GHz PC/104 processor card. Computed data are transmitted for visualization and further treatment using an Ethernet link. Main characteristics of the radar are described in Table 1.

Map construction: the R-SLAM algorithm

In the robotics domain, the simultaneous localization and mapping (SLAM) problem asks if it is possible for a mobile robot to be placed at an unknown location in an unknown environment and for the robot to build a consistent map of this environment, while simultaneously determining its location within this map. This problem has been intensively treated, particularly for autonomous navigation problems. In our application, we have been interested by the SLAM methodology to build a representation of the environment, using a succession of panoramic radar images. One interest of this approach is that this construction is "sensor-based" and avoids the use of an absolute positioning system like GPS. The majority of SLAM formulations are based on the ability of landmarks locations in the environment, assuming that these locations are time invariant. Whatever the sensor used, the most difficult problem for this approach is beacon validation and data association [START_REF] Ho | Detecting Loop Closure with Scene Sequences[END_REF]. It means that a major difficulty is (i) to confirm that the detected beacon is a valid beacon, and (ii) to associate it with a known beacon in the computed map. The data association process is an essential problem in SLAM algorithm. It has been addressed by looking for stable beacons in the environment, or using a combination of sensors (multisensorial approach) [START_REF] Zhang | An Efficient Data Association Approach to Simultaneous Localization and Map Building[END_REF]). But non-ambiguous target detection and identification based on radar signal interpretation are still open problems [START_REF] Chen | Recent advances in simultaneous localization and map-building using computer vision[END_REF]. This is particularly true in outdoor and natural environment, where landmarks which are detected by the radar cannot be modeled as simple geometric primitives. Power/range spectra contain a large amount of information mixed with various sources of noise, so detection and tracking of particular landmarks are not simple problems [START_REF] Jose | Millimetre Wave RADAR Spectra Simulation and Interpretation for Outdoor SLAM[END_REF]). The detection of landmarks in power/range spectra can be treated through constant false-alarm rate (CFAR) techniques that estimate the level of noise and the setting of a threshold [START_REF] Rohling | Radar CFAR thresholding in clutter and multiple target situations[END_REF]. But in very poor structured environment, the clutter, which is traditionally considered as noise in radar applications, must be taken into account because it can be the only available information (clutter traditionally refers to natural echoes such as ground, sea, precipitation -rain, snow or hail-, sand storms, animalsespecially birds-or insects, etc). In such situations, a detection approach based on threshold selection will fail. It is for that reason we developed a particular R-SLAM algorithm which does not treat the landmarks individually, but considers the entire radar scan. A radar image is constituted by a scan over 360° around the vehic le and is considered as an occupancy grid (matrix of cells), which divides the area of interest in regular cells. Each cell is related to a small area (20 cm × 20 cm), the occupation level of which is given by the power of the reflected signal. The global grid is updated for each 360° scan, using a data merging process. The first scan is directed at 0° referring to the initial direction of the vehicle. The absolute orientation in reference to magnetic north may be obtained by the addition of a compass. The R-SLAM process includes two main steps: -map matching between the current radar scan and the constructed map, in order to estimate the inter-scan displacement, -global map updating to integrate the current radar scan through a data merging process.

Map Matching

Classical feature extraction and matching imply to recognize the same feature repeatedly during the map construction. Then it is necessary to detect persistent and reliable features in the scan. But because of limited angle and distance resolutions, radar scans can not offer the same resolution than optic images or laser scans. Another problem is that target's reflected power depends on the orientation of the target, so the reflected signal will be different according to the angle of observation: such a property is particularly disadvantageous in performing long-term tracking of individual landmarks.

Considering the particularities of radar data, the best adapted method here is scan-matching, which performs a 3D cross correlation (two translations and one rotation), between the current image and the constructed map. The maximum of the correlation matrix gives the estimate of the radar displacement and rotation between two scans.

Global Map Updating

Once the inter-scan displacement and rotation estimated, the observation from the current radar scan has to be integrated in the global radar map.

The radar data interpretation requires a representation tool that facilitates the study of sensor-scene interaction. In that sense occupancy grids are powerful and flexible representation tools that capture the relevant scene aspects for radar-data interpretation [START_REF] Martin | Robot Evidence Grids[END_REF]. Occupancy grid represents environment by fine-grained, metric grids of variables that reflect its occupancy. Each cell of the occupancy grid contains a number that represents accumulated evidence of occupancy or emptiness in the space. In the R-SLAM process, the number contained in the cell is proportional to the microwave power of the reflected signal. This reflected power is proportional to the radar cross section (RCS) of the target. For most common types of targets such as buildings, vegetation and terrain, the radar cross section does not necessarily bear a simple relationship to the physical area, except that the larger the target size, the larger the cross section is likely to be. As the radar moves, a data merging process is applied in order to estimate the value of the cell according to the successive values of the reflected radar signal. This data merging process faces an ambivalent situation:

-the necessity to take into account the new information from the current radar image, -the necessity to preserve the information previously stored in the global map. The principle retained is to consider that the probability for a target to appear in the map is proportional to the number of times it has been observed at the same place.

Results

Two experiments are presented. Radar data are acquired and stored in real time during these experiments. The R-SLAM algorithm has been implemented in Matlab language on a Quadri-Core 2.5 GHz Pentium PC with an off-line process. It takes 3 seconds to compute a radar image and to update the radar map. Simultaneously with radar data acquisitions, GPS data are recorded in order to have a localization reference.

For the first experiment, the radar is placed on a Cemagref research vehicle, 3.3 meters above the soil. The vehicle moves on a country road in the middle of a natural and low structured extended environment near Clermont-Ferrand, France (latitude 45°44'41" N, longitude 3°07'04" E, elevation 465 m). The aerial image (from Google Earth; date: August 31st, 2006) is presented in Fig. 3(a), and the GPS trajectory is superimposed with yellow dots. The covered distance is about 1374 m, with a mean velocity of 7.9 km/h (maximum velocity: 16.6 km/h). A total of 610 radar images was stored and treated to work out the global radar map and trajectory (red dots) presented in Fig. 3(b). Few strong cooperative targets such as buildings are present within the area (lower-left part of the trajectory). Most of the elements such as trees, soil surfaces, and vegetation are slightly cooperative targets.

The second experimentation has been done on the Aydat Lake near Clermont-Ferrand, France (latitude 45°39'51" N, longitude 2°59'09" E, elevation 840 m). The radar is positioned at the front of an experimental boat from the Cemagref Institute. The radar antenna is 1.5 m above the water surface. An aerial view of the lake (from Google Earth; date: April 24th, 2004) is presented in Fig. 4(a). One can observe that a part of the lake is hidden by clouds. The covered distance is about 4 km, with a mean velocity of 1.2 m/s. Yellow dots indicate the successive GPS positions. A total of 3176 panoramic radar images was stored and processed to work out the global map and trajectory (red dots) presented in Fig. 4(b).

Conclusion

For outdoor environment perception, microwave radar is a sensor which presents a considerable potential, overcoming the limitations of vision-based sensors. In order to construct consistent maps of the environment, the R-SLAM algorithm presented in this paper exploits all the richness of radar signal, avoiding the detection of particular landmarks. The algorithm builds simultaneously a map of the environment and a localization of the vehicle within the map, without using a global positioning system. The range of radar, combined with the R-SLAM algorithm performances, makes it possible to deal with extended and low structured environments.

Figure 1 .

 1 Figure 1. Example of radar power spectrum. The frequency peaks indicate the presence of two targets (A and B) in front of the antenna.

Figure 2 .

 2 Figure 2. General view of the K2Pi radar.

  Figure 3. R-SLAM result in natural and outdoor environment. The total length of the loop is 1374 meters, the mean velocity is 7.9 km/h. (a) Aerial image from Google Earth. Yellow dots indicate GPS trajectory. (b) Global radar map and estimated trajectory (red dots). A total of 610 radar images was treated to work out the map and trajectory