J J Verbeek
email: jverbeek@science.uva.nlvlassis@science.uva.nlkrose@science.uva.nl

N Vlassis

B J A Kröse

Non-linear Feature Extraction by the Coordination of Mixture Models

Keywords: Canonical Correlation Analysis, Principal Component Analysis, Self-organizing Maps

We present a method for non-linear data projection that offers non-linear versions of Principal Component Analysis and Canonical Correlation Analysis. The data is accessed through a probabilistic mixture model only, therefore any mixture model for any type of data can be plugged in. Gaussian mixtures are one example, but mixtures of Bernoulli's to model discrete data might be used as well. The algorithm minimizes an objective function that exhibits one global optimum that can be found by finding the eigenvectors of some matrix. Experimental results on toy data and real data are provided.

Introduction

Much research has been done in the field of unsupervised non-linear feature extraction to obtain a low dimensional description of high dimensional data that preserves important properties of the data. 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. A typical example of such data are images of an object under different conditions (e.g. pose and lighting) obtained with a high resolution camera. A low dimensional representation may be desirable for different reasons like compression for storage and communication, for visualization of high dimensional data, and as a preprocessing step for further data analysis or prediction tasks.

An example is given in Fig. 1, we have data in IR 3 which lies on a two dimensional manifold (top plot). 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 'latent coordinates', i.e. coordinates on the manifold (second plot).

Here, we consider a method to integrate several local feature extractors into a single global representation. The idea is to learn a mixture model on the data where each mixture component acts as a local feature extractor. A mixture of Probabilistic Principal Component Analyzers [START_REF] Tipping | Mixtures of probabilistic principal component analysers[END_REF] is a typical example. After this mixture has been learned, the local feature extractors are 'coordinated' by finding for each local feature extractor a suitable linear map into a single 'global' low-dimensional coordinate system. The collective of the local feature extractors together with the linear maps into the global latent space provides a global non-linear projection from the data space into the latent space.

This differs from other non-linear feature extraction methods like: Generative Topographic Map (GTM) [START_REF] Bishop | GTM: The generative topographic mapping[END_REF], Locally Linear GTM [START_REF] Verbeek | Locally linear generative topographic mapping[END_REF], Kohonen's Self-Organizing Map [START_REF] Kohonen | Self-Organizing Maps[END_REF] and a Variational Expectation Maximization algorithm similar to SOM [START_REF] Verbeek | Self-organization by optimizing free-energy[END_REF]. These algorithms are actually quite opposite in nature to the method presented here: they start with a given configuration of mixture components in the latent space and adapt the mixture configuration in the data space to optimize some objective function. Here we start with a given and fixed mixture in the data space and find an optimal configuration in the latent space. The advantage of the latter method is that it is not prone to local optima in the objective function. 1 Moreover, the solution for a d-dimensional latent space can be found by finding the eigenvec-tors of a matrix corresponding to the 2 nd smallest up to the (d + 1) st smallest eigenvalues. The size of this matrix is given by the total number of locally extracted features plus the number of mixture components, i.e. if all k components extract d features then the size is k(d + 1).

The work presented in this paper can be categorized together with other recent algorithms for unsupervised non-linear feature extraction like IsoMap [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF], Local Linear Embedding [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF], Kernel PCA [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF], Kernel CCA [START_REF] Bach | Kernel independent component analysis[END_REF] and Laplacian Eigenmaps [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF]. All these algorithms perform nonlinear feature extraction by minimizing an objective function that exhibits a limited amount of stable points that can be characterized as eigenvectors of some matrix. These algorithms are genrally simple to implement, since it is only needed to construct some matrix, then the eigenvectors are found by standard methods.

In this paper we build upon recent work of Brand [START_REF] Brand | Charting a manifold[END_REF]. In the next section we review the work of Brand, we derive the objective function by a slight modification of the free-energy used in [START_REF] Roweis | Global coordination of local linear models[END_REF] and we discuss how we can find its stable points. In Section 3 we present a method to obtain better results by using several mixtures for the data in parallel. Secondly, we show in Section 4 that finding the locally valid linear maps between the data space and the latent space can be regarded as a weighted form of Canonical Correlation Analysis (CCA), and we show how we can generalize the work of Brand to perform nonlinear CCA. A general discussion and conclusions can be found in Section 5.

Linear Embedding of Local Feature Extractors

In this section we describe the work of Brand [START_REF] Brand | Charting a manifold[END_REF], in a notation that allows us to draw links with CCA in section 4.

Consider a given data set X = {x 1 , . . . , x N } and a collection of k local feature extractors of the data f s (x) for s = 1, . . . k. Thus, f s (x) is a vector containing the features extracted from x by feature extractor s. Each feature extractor gives zero or more (linear) features of the data. With each feature extractor, we associate a mixture component in a probabilistic k-component mixture model. We write the density on high dimensional data items x as:

p(x) = k s=1 p(s)p(x|s), (1)
where s is a variable ranging over the mixture components and p(s) is the a-priori probability on component s (sometimes also termed 'mixing weight' of component s).

Next, we consider the relationship between the given representation of the data (denoted with x) and the representation of the data in the latent space, which we have to find. Throughout, we will use g to denote latent coordinates for data items. The g is for 'Global' latent coordinate as opposed to locally extracted features, which we will denote by z. For the unobserved latent coordinate g corresponding to a data point x and conditioned on s, we assume the density:

p(g|x, s) = N (g; κ s + A s f s (x), σ 2 I) (2)
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 mean κ 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 marginal (over s) distribution on latent coordinates given data:

p(g|x) = s p(s, g|x) = s p(s|x)p(g|x, s). (3)
Given a fixed set of local feature extractors and a corresponding mixture model, we are interested in finding linear maps L s that give rise to 'nice' projections of the data in the latent space. By 'nice', we mean that components that have a high posterior probability to have generated x (i.e. p(s|x) is large) should give similar projections of x in the latent space. So we would like the p(g|x, s) to look similar for the components with large posterior. Since the marginal p(g|x) is a Mixture of Gaussians, we can measure the similarity between the predictions by looking how much the mixture resembles a single Gaussian. If the predictions are all the same, the mixture is a single Gaussian.

We are now ready to formally define the objective function Φ({L 1 , . . . , L k }), which sums the data log-likelihood and a Kullback-Leibler divergence D which measures how uni modal the distribution p(g|x) = s p(s|x)p(g|x, s) is.

Φ = N n=1 log p(x n ) - n D(Q n (g, s) p(g, s|x n )), (4)
where Q n distributes independently over g and s:

Q n (g, s) = q ns Q n (g) = q ns N (g; g n , Σ n ). ( 5
)
This objective function is similar to the one used in [START_REF] Roweis | Global coordination of local linear models[END_REF][START_REF] Verbeek | Coordinating Principal Component Analyzers[END_REF]. The difference here is that the covariance matrix of p(g|x, s) does not depend on the linear maps L s . If we fix the data space model, turning the data log-likelihood into a constant, and set q ns = p(s|x n ), then the objective sums for each data point x n the Kullback-Leibler divergence between Q n (g) and p(g|x, s), weighted by the posterior p(s|x n ):

Φ = n,s q ns D(Q n (g) p(g|x n , s)). (6)
In order to minimize the objective Φ with respect to g n and Σ n , it is easy to derive that we get:

g n = s q ns g ns and Σ n = σ 2 I, (7)
where I denotes the indentity 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 ≥ 0, (8)
with g ns = L s z ns as before. Note that we can rewrite the objective as:

Φ = 1 2 n,s,t q ns q nt g nt -g ns 2 ≥ 0. ( 9
)
The objective is a quadratic function of the linear maps L s [START_REF] Brand | Charting a manifold[END_REF]. At the expense of 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 ], (10)
U = [u 1 . . . u N ] , (11)
L = [L 1 . . . L k ] . (12)
Note that from ( 7), [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF] and [START_REF] Roweis | Global coordination of local linear models[END_REF] we have: g n = (u n L) . The expected projection coordinates can thus be computed as:

G = [g 1 . . . g N ] = UL. ( 13
)
We define the block-diagonal matrix D as:

D =     D 1 0 . . . 0 D k     (14)
where D s = n q ns z ns z ns . The objective then can be written as:

Φ = Tr{L (D -U U)L}. ( 15
)
The objective function is invariant for translation and rotation of the global latent space, as can be easily seen from [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF]. Furthermore, re-scaling the latent space changes the objective monotonically.

To make solutions unique with respect to translation, rotation and scaling, we put two constraints:

(trans.) : ḡ = 1 N n g n = 0, (16)
(rot. + sc.) :

Σ g = n (g n -ḡ)(g n -ḡ) (17) = L U (I -1/N )UL = I,
where we used 1 to denote the matrix with ones everywhere. Solutions are found by differentiation using Lagrange multipliers, and the columns v of L are found to be generalized eigenvectors:

(D -U U)v = λU (I -1/n)Uv (18)
Since the solutions are zero mean, it follows that 1Uv = 0, and hence that:

(D -U U)v = λU Uv (19) ⇔ Dv = (λ + 1)U Uv. ( 20
)
The value of the objective function is given by:

Φ = i v i (D -U U)v i (21) = i λ i v i U Uv i (22) = (N -1) i λ i var(g i ), (23)
where var(g i ) denotes the variance in the i th coordinate of the projected points. So rescaling to get unit variance makes the objective equal to (N -1) i λ i . The smallest eigenvalue is always zero, corresponding to projections that collapses all projections in the point κ. This projection has for s = 1, . . . , k : A s = 0 and κ s = κ. This is embedding tells us nothing about the data, therefore we need the eigenvectors corresponding to the second up to the (d + 1) st smallest eigenvalues to obtain the best embedding in d dimensions.

Note that there is no problem if different feature extractors extract a different number of features.

In Fig. 1 we give an illustration of the above. The top two images show the original data presented to the algorithm and the 2-dimensional representation found by the algorithm. The blue and yellow sticks indicate the local feature extractors which are centered on the red dots. The two lower scatter plots compare the coordinates on the data generating manifold with the discovered latent coordinates. For both coordinates there is high correlation with the true latent coordinates. 3 Incorporating more points in the objective function

In this section we discuss how using multiple mixtures for the data can help obtain more stable results. The next section then discusses how this method can be exploited to define a non-linear form of Canonical Correlation Analysis.

Sometimes the method of the previous section collapses several linear feature extractors in some part of the underlying data manifold onto a point or a line. Whereas the Automatic Alignment (AA) method [START_REF] Teh | Automatic alignment of local representations[END_REF] was less sensitive to such phenomena.

AA is based on the Local Linear Embedding algorithm [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] and finds linear maps to integrate a collection of locally valid feature extractors in a manner similar to what we described in the previous section. In AA the objective function is:

Φ AA = m,n W nm g n -g m 2 , (24)
again we have g n = s q ns g s n . The weights, collected in the weight matrix W, are found by first determining for every data point x n its nearest neighbors nn(n) in the data space and use the W that minimizes:

n x n - m∈nn(n) W nm x m 2 , (25)
under the constraint that m W nm = 1 and only nearest neighbors can have non-zero weights, thus: ∀m / ∈ nn(n) : W nm = 0. The objective function [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF] of the previous section is in fact only based on points for which the posteriors q ns assigns non-negligible mass on at least two mixture components. This is seen directly from ( 9), since if one of the q ns takes all mass, all products of the posteriors are zero except for one term where the squared distance is zero. In other words, the objective only focuses on data that is on the overlap of at least two mixture components. So, only points for which the entropy of the posterior is reasonably high contribute to the objective function. The other data, for which the entropy in the posterior is low, is thus completely neglected in the objective function and thus in the coordination of the local feature extractors.

In practice, often the bulk of the data has a low entropy posterior and hardly plays a role in the objective function. To make more use of the available data we should have more entropy in the posteriors. One way to do that is to smoothen the posteriors by taking the closest posterior in KL-divergence sense to the old posterior that has a certain entropy, yielding: q ns ∝ q α ns . Another possibility is to smoothen such that the average entropy has a certain value. However, to determine the 'right' amount of entropy in the posteriors seems a rather problematic question.

The problem can also be circumvented by using several mixtures, say two, in parallel. The mixtures are trained independently and the final posterior for a specific mixture component from one of the two mixtures is defined as half of the posterior within the mixture from which this component originates. So, the sum of the posteriors on the components originating from the first mixture equals 1/2 and similarly for the second mixture. In this manner we always have an entropy of at least one bit in the posteriors, and hence all points make a significant contribution to the objective function.

Note that we are using just a heuristic here, more principled ways to learn a mixture with high entropy in the posteriors could be divised but do not seem essential. Only in the case that the two independently learned mixtures are the same we have not gained anything. This, however, is not likely to be the case for mixture models fitted by EM algorithms that find local optima of the loglikelihood function. We initialize the EM mixture learning algorithm at different values and in this manner we (probably) find mixtures that are different local optima of the data log-likelihood function. We can think of this method as covering the data manifold not with one set of tiles, but using two coverings of the manifold. In this way every tile is overlapping with several other tiles on its complete extent.

We found in experiments that using two mixtures instead of one gives great improvement in the obtained results. Moreover, we prevent the computation of the weight matrix W needed by AA. The time needed to compute the weight matrix is in principle quadratic in the number of data points and can be problematic in non-Euclidean spaces. By using several mixtures in parallel, we avoid the quadratic cost and have a method that is based purely on mixture models and that can be readily applied to any type of mixture models.

In Fig. 2 we give some examples of data embeddings obtained on a data set of 1965 pictures of a face. 2 First the images were projected from the original 28 × 20 pixel coordinates to 10 dimensions with PCA. The projection from 10 to 2 dimensions was done by fitting a 10 component mixture of 2-dimensional probabilistic PCA [START_REF] Tipping | Mixtures of probabilistic principal component analysers[END_REF] and aligning them with the method described in the previous section. We show the latent co-2 These images can be obtained from the web page of Sam Roweis at http://www.cs.toronto.edu/~roweis. ordinates assigned to the pictures as dots, some pictures are shown next to their latent coordinate.

The two bottom figures, where we used two parallel mixtures, show a nice separation of the 'smiling' faces and the 'angry' faces. Furthermore, the variations in gaze direction are nicely identified. The two top images, where we used single mixtures with respectively ten and twenty mixture components do not show the separation between the two 'moods'.

Non-linear Canonical Correlation Analysis

In Canonical Correlation Analysis (CCA) [START_REF] Hotelling | Relation betweeen two sets of variates[END_REF][START_REF] Mardia | Multivariate Analysis[END_REF] 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 with minimizing:

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

[ n x n x n ]a + b[ n y n y n ]b = 1. Several generalizations of CCA are possible.
Generalizing the above such that the sets do not need to be zero mean and allowing a translation as well gives the problem of minimizing:

E = 1 2 n [(ax n + t a ) -(by n + t b )] 2 , (27)
again under the constraints that the sum of variances of the projections of X and Y are equal to some fixed constant. We can also generalize by mapping to IR d instead of the real line, and then requiring the covariance matrix of the projections to be identity. 3 CCA can also be readily extended to take into account more than two point sets, c.f. [START_REF] Bach | Kernel independent component analysis[END_REF].

In the generalized CCA with multiple pointsets and allowing translations and mapping to IR d , we minimize the squared distance between all pairs of projections under the constraints that each set of projected points has unit variance. We denote the projection of the n-th point in the sth point-set as g ns . We thus minimize the error 3 If the linear maps are constrained to be orthonormal matrices (rotations plus reflection) then this is known as the Procrustes rotation [START_REF] Cox | Multidimensional Scaling[END_REF]. function:

Φ CCA = 1 2k n k s=1 k t=1 g ns -g nt 2 (28) = n s g ns -g n 2 , ( 29
)
where g n = 1 k s g ns . The constraint is that:

s n g ns 2 = 1.
Note that the objective Φ in equation ( 8) coincides with Φ CCA in (29) if we set q ns = 1/k. The different constraints imposed upon the optimization by CCA and our objective of the previous sections turn out to be equivalent, since the stable points under both constraints coincide. We can regard the features extracted by each local feature extractor as a point set. Hence, we can interpret the method of the Section 2 as a weighted form of CCA that allows translation on multiple point sets.

Viewing the coordination of local feature extractors as a form of CCA suggests using the coordination technique for non-linear CCA. This is achieved quite easily, without modifying the objective function [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF]. We consider different point sets, each having a mixture of locally valid linear projections into the 'global' (in the sense of shared by all mixture components and point sets) latent space. 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. As with the use of parallel mixtures, the weight q ns associated with the s th projection of observation n, is given by the posterior of the mixture component in the corresponding mixture (e.g. the posterior p(s|x n ) if projection s is due to a component in the mixture for X and p(s|y n ) if s refers to a component in the mixture fitted to Y) divided over the number of point sets. It is now clear that the use of parallel mixtures in the previous section was merely a form of non-linear CCA where both point sets are the same.

We again define g n = s q ns g ns , where s now ranges over all mixture components in all mixtures fitted on the C different point sets. We use q c nr to denote the posterior on component r for observation n in point set c. where the q c nr are rescaled with a factor C such that all q c nr for one point set c sum to one. Furthermore, we define g c = r q c nr g c nr as the average projection due to point set c. The objective can be rewritten in different forms:

Φ = 1 2 n s,t q ns q nt g ns -g nt 2 (30) = n s q ns g ns -g n 2 (31) = c 1 C n,s q c ns g n -g c ns 2 (32) = c 1 C n g n -g c n 2 + c 1 C n s q c ns g c n -g c ns 2 . ( 33
)
Observe how in (33) the objective sums both between point set consistency of the projections (first summand) and within point set consistency of the projections (second summand).

As an illustrative toy application of the nonlinear CCA we took two point-sets in IR 2 of 600 points each. The first point-set was generated on an S-shaped curve the second point set was generated along a circle segment. To both point sets we added Gaussian noise. We learned a k = 10 component mixture model on both point-sets, illustrated in the top figures of Fig. 3. In the bottom figure the discovered latent space representation (vertical axis) is plotted against the coordinate on the generating curve (horizontal axis). The relationship is monotonic and roughly linear.

Discussion and Conclusions

We have shown how we can globally maximize a free-energy objective similar to that of [START_REF] Roweis | Global coordination of local linear models[END_REF][START_REF] Verbeek | Coordinating Principal Component Analyzers[END_REF] with respect to the linear maps L s that map homogeneous coordinates of the local feature extractors into a global latent space. This was achieved by making the variance of p(g|s, x) independent of the linear maps. The linear maps are found by computing the d+1 smallest eigenvectors of a matrix. The edge size of the eigenproblem is given by k plus the total number of extracted features. Typically, if all feature extractors give the same number of features, this will be k(d+1). However, it is not necessary to let every mixture component extract an equal number of features.

The way in which we find the latent representation of the data is very similar to the k-segments algorithm [START_REF] Verbeek | A k-segments algorithm for finding principal curves[END_REF]. In the latter we first fitted a mixture of PCA's and used a (only locally optimal) combinatorial optimization algorithm to discover how the different PCA's (segments) could be aligned. Also here, first a mixture model is fitted to the data and then an optimization problem is solved to coordinate the mixture components. However, in the current work the optimization of the alignment is globally optimal and can be applied to PCA's of any dimension and latent spaces of any dimension.

As compared to [START_REF] Teh | Automatic alignment of local representations[END_REF] we do not have to compute neighbors anymore in the data space, which is costly and moreover can be hard when processing data with mixed discrete and continuous features. We only need a mixture model to produce posteriors here. In [START_REF] Brand | Charting a manifold[END_REF] the same objective and solution are used as here. We showed how using multiple mixtures in parallel can boost performance significantly and how the same technique can be used for non-linear CCA. We illustrated the viability of the presented algorithms on synthetic and real data.

In future work, we want to exploit the link with CCA further by investigating its applicability in classification and regression problems.

Figure 1 :

 1 Figure 1: From top to bottom. (1) Data in IR 3 with charts indicated by the blue and yellow sticks. (2) Data representation in IR 2 . (3) First original coordinate on manifold against first found latent coordinate. (4) Second original coordinate on manifold against second found latent coordinate.

Figure 2 :

 2 Figure 2: From top to bottom: (1) One mixture of 10 components. (2) One mixture of 20 components. (3) Two parallel mixtures of 10 components each.

Figure 3 :

 3 Figure 3: (Top) Two point-sets and the fitted mixture models. Local charts are indicated by the blue sticks. (Bottom) Scatter plot of the found latent coordinates (vertical) and coordinate along the generating curve (horizontal).

Of course, fitting the mixture model on the data might be susceptible to local optima.

Acknowledgment

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


Josep M Porta
email: porta@science.uva.nl

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

Ben Kröse
email: krose¡@science.uva.nl

Enhancing Appearance-based Robot Localization Using Sparse Disparity Maps

In this paper, we enhance appearance-based robot localization by using disparity maps. Disparity maps provide the same type of information as range based sensors (distance to objects) and thus, they are likely to be less sensitive to changes of illumination than plain images, that are the source of information generally used in appearancebased localization. The main drawback of disparity maps is that they can include very noisy depth values: points for which the algorithms can not determine reliable depth information. These noisy values have to be discarded resulting in missing values. The presence of missing values makes Principal Component Analysis (the standard method used to compress images in the appearance-based framework) unfeasible. We describe a novel Expectation-Maximization algorithm to determine the principal components of a data set including missing values and we apply it to disparity maps. The results we present show that disparity maps are a valid alternative to increase the robustness of appearancebased localization.

I. INTRODUCTION

In appearance-based localization [START_REF] Kröse | A Probabilistic Model for Appearancebased Robot Localization[END_REF], the environment is not represented explicitly, but implicitly as a database of images (or features derived from them) collected at known positions in a training phase. For localization, the features of the image observed by the robot are compared with the features stored in the database.

A problem is that images largely change due to variations on illumination and so do the corresponding features. Consequently, illumination changes can largely degrade the performance of the appearance-based localization techniques. A possible solution to this problem is to complement intensity information with information provided by range sensors such as sonar or laser scanners.

In this paper, we describe a method which uses the depth information from a stereo vision system to enhance the appearance-based localization. Depth profiles, or disparity maps, are computed by matching points on two images. If the depth is only computed for few of the points in the images we obtain a non-dense disparity map that do not provide much more information that that provided by sonars or laser scanners. If we compute the depth for all pixels then we get a dense disparity map [START_REF] Alvarez | Dense Disparity Map Estimation Respecting Image Discontinuities: a PDE and Scale-Space Based Approach[END_REF]. The advantage of this kind of maps is that, in contrast to standard laser range finders, which provide a 1-D depth map, a more informative 2-D depth map is determined.

Additionally, dense disparity maps could be directly used for appearance-based localization as if they were plain images. The problem of dense disparity maps is that they can include very noisy depth values (especially in nontextured areas in the scene). If these non-reliable depth values are taken apart, we get a sparse disparity map, that is, a disparity maps that include some missing values. The issue of missing values particularly plays a role in the way we extract features for the appearance model. For the intensity images, we use features extracted from a Principal Component Analysis (PCA) on the set of collected images [START_REF] Kröse | A Probabilistic Model for Appearancebased Robot Localization[END_REF]. The numerically most accurate way to compute the features from the image set is by taking the Singular Value Decomposition (SVD) of the image matrix. However, because of the missing values, this technique can not be applied to sparse disparity maps.

In this paper, we propose the use of a Expectation-Maximization (EM) algorithm to determine the principal components for a set of disparity maps including missing values. This algorithm is similar to that described in [START_REF] Roweis | EM Algorithms for PCA and SPCA[END_REF], but we introduce a simpler way to deal with the missing values that makes the whole process more efficient. Thanks to this EM algorithm sparse disparity maps can be used in appearance-based localization.

This paper is organized as follows. First, we introduce the general framework for appearance-based localization. Next, we briefly describe how disparity maps are defined from intensity images. In Sections IV to VI, we introduce the EM algorithm to extract features from the disparity maps and we describe how these features can be used in our appearance-based localization framework in addition to features derived from intensity images. Next, we present results obtained in a real environment and, finally, we conclude summarizing our work and extracting some conclusions out of it.

II. APPEARANCE-BASED LOCALIZATION

Appearance-based localization departs from a set of training images Z ¢ ¤£ z 1 ¥ §¦ §¦ ¨¦ §¥ z N © taken at known poses (position/orientation) X ¢ £ x 1 ¥ §¦ ¨¦ §¦ ¥ x N © . Images are linearly compressed to get a set of features for each image Y ¢ £ y 1 ¥ §¦ ¨¦ §¦ ¥ y N © , where

y i ¢ Cz i ¦ (1)
The projection matrix C expands the principal components of the data set Z and it is computed off-line applying PCA. The rows of C are the eigenvectors corresponding to the d largest eigenvalues.

In the on-line execution, we aim at estimating the probability on the robot's pose at time t, p£ x t © . This is usually done assuming that the environment is Markovian and deriving p£ x t © from p£ x t 1 © , from a known action model, p£ x t u t 1 ¥ x t 1 © (with u t 1 the action executed by the robot), and from a sensor model, p£ y t x t © . In the appearance-based framework, the sensor model is defined using the training set and the current observation, y t . For instance, Vlassis et al. in [START_REF] Vlassis | Auxiliary Particle Filter Robot Localization from High-Dimensional Sensor Observations[END_REF] propose a nearest-neighbor formulation for the sensor model

p£ y x © ¢ J ∑ j 1 λ j φ £ x x j © ¥ (2)
with x j the J nearest neighbors (i.e., the subset of training points obtained at positions x j and with an associated set of features y j more similar to the features of the current observation y t ), λ j a set of weights that favor closer nearest neighbors, and φ a Gaussian. The closer the training points used to define the sensor model to the actual position of the robot, the smaller the error in localization. Changes in illumination result in wrong matches between the features of the current image and those in the training set and, thus, in a wrong sensor model and in a wrong update of the robot's position. To minimize this problem as much as possible, we propose to make use of the information provided by disparity maps.

III. DISPARITY MAPS

Disparity maps are determined matching points in images taken by a pair of calibrated cameras mounted on the robot. Given a single image, the three-dimensional location of any visible object point Q must lie on the line that passes through the center of projection of the camera c and the image of the object point p. The determination of the intersection of two such lines generated from two independent images is called triangulation and provides the 3-D position of Q w.r.t the cameras.

To apply triangulation, for each pixel p in one of the images we have to search for a correspondent point in the other image. The epipolar constraint reduces this search to a single line called the epipolar line [START_REF] Faugeras | Three-Dimensional Computer Vision: A Geometric Viewpoint[END_REF]. Usually, the correspondence is done by comparing areas around pixel p with areas around each candidate pixel p . The most similar pixels p and p are assumed to correspond to different projections of the same point Q in the scene. If the images planes for the two cameras are co-planar, the distance r from the scene point Q to the cameras can be computed as

r ¢ b f d d ¥ PSfrag replacements Q c c d d p p f f b r Fig. 1. Elements in disparity computation.
where b is the baseline (distance between the two viewpoints), f is the focal length of the cameras, d is the horizontal distance from the projected point p to the center of one of the images, and d is the same for the other image (see Figure 1). The difference d d is called disparity.

With constant baseline b and focal length f , the disparity is inversely proportional to the depth of point Q.

The stereo algorithm we use [START_REF] Konolige | Small Vision System: Hardware and Implementation[END_REF] applies many filters in the process to determine the disparity map both to speed up the process and to ensure the quality of the results. For instance, if the area around pixel p is not textured enough it would be very difficult to find a single corresponding point p (we are more likely to end up with many points p with almost the same probability of being the corresponding point of p). For this reason, pixels on low textured areas are not even considered in the matching process. The result of this and other filtering processes to avoid too noisy disparity values is to produce a sparse disparity map: a disparity map where many pixels do not have a disparity value (see Figure 2). This makes the use of standard PCA to determine the projection matrix unfeasible and we have to use more elaborated techniques as the EM algorithm introduced in the following sections.

IV. THE EM ALGORITHM

The EM algorithm [START_REF] Dempster | Maximum-Likelihood for Incomplete Data via the EM Algorithm[END_REF], [START_REF] Minka | Expected-Maximization as Lower Bound Maximization[END_REF] aims at determining the maximum likelihood model for a given set of data Z. Not all possible models are considered, but only those inside a family of models defined by a set of parameters, collectively denoted as θ (in our application, these parameters are a projection matrix C and a parameter σ that indicates the noise in the data out of the principal component subspace). The optimal model is determined by maximizing the likelihood or, equivalently, the log-likelihood: 1 Φ£ θ © ¢ log p£ Z; θ © ¦ In many cases, Φ can be defined using an auxiliary set of unobserved or hidden variables h. In this case, we have 1 The log is used since it makes the resulting expressions simpler. Φ might be maximized using a gradient ascent process.

Φ£ θ © ¢ log h p£ Z ¥ h; θ © ¦
However, the EM algorithm provides a simpler to implement and step-size free algorithm using iterative lowerbound maximization.

For a given set of parameters θ , the EM algorithm first finds a lower bound Ψ of Φ, possibly such that Ψ touches Φ at point θ . After that, the lower bound Ψ is maximized for the parameters θ . The definition of Ψ for a fixed set of parameters θ is the E step of the algorithm and the maximization of the resulting Ψ is called the M step. The sequential iteration of E and M steps from an arbitrary initial set of parameters θ is guaranteed to find a (local) maximum of Φ.

The lower-bound Ψ used in the E step, can be defined as

Ψ£ θ © ¢ Φ£ θ © KL£ q£ h © p£ h Z; θ © §© "! Φ£ θ © ¥
where KL denotes the (non-negative) Kullback-Leibler divergence between two probability distributions. To make Ψ£ θ © ¢ Φ£ θ © for the current parameters θ , the KL should be zero. The KL divergence is zero iff the two compared probability distributions are equal. Thus, in the E step we set q£ h © ¢ p£ h Z; θ © .

In the M step, we have to maximize Ψ. We can rewrite Ψ as

Ψ£ θ © ¢ Φ£ θ © # q£ h © q£ h © p£ h Z; θ © ¢ ¢ $ E q% h& "' log q£ h © p£ h Z; θ © )( 10 log p£ Z; θ © ¢ ¢ E q% h& 2' log p£ Z ¥ h; θ © q£ h © ( ¢ ¢ h q£ h © log p£ Z ¥ h; θ © q£ h © logq£ h © ¦
The relevant term of Ψ to be maximized is

Ψ M ¢ 3 h q£ h © log p£ Z ¥ h; θ ©
since the other term of Ψ does not depend on θ .

For a more detailed presentation of the EM algorithm see [START_REF] Minka | Expected-Maximization as Lower Bound Maximization[END_REF], which we followed in the above brief description.

V. APPLYING EM TO DISPARITY MAPS

Given a D 4 N matrix Z, representing a set of disparity maps Z ¢ 65 z 1 ¥ §¦ ¨¦ §¦ ¥ z N 7 , our aim is to find a projection that maps Z onto its d-dimensional principal components space (d 8 D). For each disparity map z i , we have a subset of elements that are observed z i9 o and a set of missing values z i9 h . Thus, we want to define the probability p£ y i z i9 o © ¥ with y i the set of principal component features derived from the observed values in disparity map z i . We do so by first defining p£ y © as a Gaussian with zero mean and identity covariance, and p£ z y © ¢ A@ B£ z;CC y ¥ σ 2 I © . The projection density p£ y i z i9 o © then follows a Gaussian distribution @ D£ y i ; µ y i ¥ Σ y i © , where Σ 1

y i ¢ I 0 C o CC o E σ 2 and µ y i ¢ Σ y i C o z o E σ 2 .
Here, C is the d 4 D projection matrix from the space of disparity maps to the space of principal components, CC performs the inverse mapping, and C o contains only the columns of C that correspond the observed elements in a given map z i . Observe that C o is different for each disparity map and hence so is the covariance matrix Σ y i . As mentioned, the parameters that describe our model are the projection matrix C and σ .

To find the maximum likelihood projection matrix C, we have to particularize the general E and M steps described in the previous section. This can be done paying attention only to the observed values z o , taking the projected features as the hidden variables, and not making inference about the missing pixels. However, the EM algorithm can be made more efficient by estimating y and z h at the same time.

In the E step the KL divergence to be minimized is

KL ¢ F q£ y ¥ z h © q£ y ¥ z h © p£ z h ¥ y z o ;C ¥ σ © ¦
In order to obtain an efficient E step we factor q over y and z h . Thus,

KL ¢ 3 q£ y © ' log q£ y © p£ y z o ;C ¥ σ © 0 q£ z h © log q£ z h © p£ z h y;C ¥ σ © G( ¦
Using free-form optimization, we have

q£ z h © ∝ exp q£ y © log p£ z h y;C ¥ σ © ¥ q£ y © ∝ p£ y z o ;C ¥ σ © exp q£ z h © log p£ z h y;C ¥ σ © ¦
Both of these densities are Gaussian. The covariance of the q£ z h © is given by σ 2 I and the mean for the missing values on the disparity map i is given by z i9 h ¢ CC h y i . The covariance and the means (collected as columns in Y ) of the q£ y © for the disparity maps are given respectively by

Σ y ¢ $H I 0 σ 2 CCC PI 1 ¥ Y ¢ σ 2 Σ y CZ ¥
where Z is the data matrix with missing values filled in with the new z i9 h . Observe that, the covariance matrix Σ y is common for all data points when computing the set of projections Y . Therefore, we only have to perform one matrix inversion and not one per data point.

With respect to the M step, the relevant part of the likelihood function to be optimized is

Ψ M ¢ 3 y9 z h q£ y ¥ z h © log p£ z o ¥ z h ¥ y;C ¥ σ © ¦
Analytical manipulation leads to the effective objective to be maximized:

Ψ M ¢ ND logσ 2 0 0 σ 2 Q N ∑ i 1 z i CC y i 2 0 Tr5 CC Σ y C 7 SR 0 0 σ 2 D h σ 2 old ¥
with D h the total amount of missing values in Z and σ old the previous value for σ . From the above, we get the simple updates

C ¢ ZY C UT NΣ y 0 YY C WV 1 σ 2 ¢ 1 ND Q NTr5 CC Σ y C 7 0 N ∑ i 1 z i CC y i 2 0 D h σ 2 oldR ¦
The E and M steps have to be iterated while there is a large (relative) change in the lower bound of the loglikelihood function that, after the updates, reads to

Ψ£ C ¥ σ © ¢ $ N 2 X D log σ 2 0 Tr5 Σ y 7 log Σ y Y 1 2 Tr5 YY C 7 0 1 2 D h log σ 2 old ¦
A reasonable initialization of C would be that containing d randomly selected disparity maps (with zeros in the missing values) and σ 2 equal to the initial reconstruction error (for the observed values).

Each iteration of the EM steps scales with order O£ dDN © , assuming d `N. Our approach to finding the maximum likelihood C and σ has the advantage that the guess of the missing values comes essentially for free while in other approaches [START_REF] Roweis | EM Algorithms for PCA and SPCA[END_REF] they were obtained extending the E step with a optimization process in the space of images for every image separately. This is an important saving especially when working with data sets such as space disparity map where each point z contains many missing values.

Once C and σ are determined, the set of features y corresponding to a new disparity map z are computed as

y ¢ a£ σ 2 I 0 C o CC o © 1 C o z o ¦ The matrix £ σ 2 I 0 C o CC o © 1 C o maps from (observed)
disparity values to feature detectors and plays the same role as matrix C in equation 1. For vanishing sigma,

£ σ 2 I 0 C o CC o © 1 C o is the pseudo inverse of C o .
The most expensive step in the computation of the vector of feature detectors for each disparity map is the product C o CC o that is O£ Dd 2 © . Because of this, the on-line determination of the features is only feasible for relative small d (and moderated D's).

A Matlab version of the EM algorithm just described can be downloaded from [START_REF] Verbeek | Probabilistic PCA with Missing Values[END_REF].

The just described algorithm can be applied to any set of data with missing values. Thus, we can apply it to disparity maps but we can also apply it to a concatenation of the intensity image and the disparity map obtained at a given position. The advantage of concatenating intensity and disparity maps is that the information from the intensity images can eventually help to deal with missing values in the disparity maps. The problem of proceeding this way is that changes in illumination would largely modify part of the input vector (that relative to intensity) and, as a consequence, we would lose the robustness to changes in illumination offered by disparity maps.

VI. SENSOR MODEL FUSION

Once we have a way to define features from disparity maps, it remains the question of how to combine the information coming from disparity with that obtained from intensity to define a unified sensor model. Two possible solutions come to mind: to combine them in a conjunctive way or in a disjunctive one.

A conjunctive-like combination can be achieved factorizing the sensor model p£ y d ¥ y i x © ¢ p£ y d x © p£ y i x © ¥ with y d the features obtained for disparity and y i those for intensity. In this way, only those training points consistent with both the current intensity image and the current disparity map are taken into account to update the robot's position. The problem of this formulation is that wrong matches for intensity (or for disparity) would result in an almost null sensor model and, thus, the position of the robot would be updated almost without sensory information. To avoid this, we propose to use a disjunctive-like model that can be implemented defining the global sensor model as linear combination of the intensity and disparity sensors models p£ y d ¥ y i x © ¢ w d p£ y d x © 0 w i p£ y i x © ¥ that, using equation 2, reads to

p£ y d ¥ y i x © ¢ w d J ∑ j 1 λ j φ £ x x j © 0 w i Jb ∑ j 1 λ j φ £ x x j © ¥
where x j and x j are the closest training points to the features of the current disparity and intensity images respectively. The weights w d and w i (with 0 ! w d ! 1, w i ¢ 1 w d ) can be used to balance the importance of the information obtained with each type of sensor. If both information sources are assumed to be equally reliable, we can set

w d ¢ w i ¢ 1 E 2.
With this expression for the sensor model, all the likely hypotheses on the robot position suggested by the current observation are taken into account. The particle filter we use to update the probability on the robot's position [START_REF] Pitt | Filtering via Simulation: Auxiliary Particle Filters[END_REF], [START_REF] Vlassis | Auxiliary Particle Filter Robot Localization from High-Dimensional Sensor Observations[END_REF] takes care of filtering the noise and, thus, of preserving the hypothesis that is more consistent over time.

VII. EXPERIMENTS AND RESULTS

To assess the contribution of using disparity maps in appearance-based localization, we collected a training set moving our robot (see Figure 3) in an area of 900 4 500 centimeters. Training points were disposed on a regular grip with a resolution of 50 cm. At each point on the grid, we took intensity images and disparity maps rotating the camera every 10c . This makes a total amount of more than 4500 training images. Each image (and each disparity We use N ¢ 100 randomly sampled images to compute the d ¢ 20 principal components of the training set. For intensity images, we used a standard PCA technique to determine the principal components and, for disparity maps, we applied the algorithm we describe in Section V.

The principal component of the disparity maps are computed in 13 iterations of our EM-based algorithm with a convergence threshold (relative change on two consecutive maximums of functions Ψ) of 0 ¦ 001. The execution of this process takes 125 seconds in a Pentium 4 at 2.0 GHz.

To test the our localization system, we moved the robot along a pre-defined path in three different illumination conditions: using tube lights (that was the light used to get the training set), using bulb lights and, finally, with natural light (obtained opening the curtains of the windows all along one of the sides of the lab). These three different illuminations ensure not only a difference in the global intensity of the light, but also a change in the distribution of light sources. At regular distances along the test path, we take an image and we compute the corresponding sensor model using the J training points with a set of features more similar to those corresponding to the just obtained image. As mentioned before, the closer the training points used to define the sensor model to the actual position of the robot, the better the sensor model and, thus, better the update of the robot position estimation.

Figure 4 shows the average (and the variance) of the error in the sensor model e ¢ min d n r n ¥ with r ¢ 3£ r x ¥ r y ¥ r φ © the pose of the robot at the test position and n ¢ e£ n x ¥ n y ¥ n φ © the pose of the points used to define the sensor model (that are different for each test position).

An error in the range H 25 ¥ 50I is quite reasonable since the distance between training points in X and Y dimensions is 50 cm.

We repeat the test in two cases: (a) using only intensity images (dashed line on Figure 4) and (b) using, additionally, disparity maps (dotted line on the figure). In the first case we use J ¢ 10 (quite small compared with the total number of images in the training set) and in the second case we use J ¢ 5 but for both intensity and disparity (so, we also get 10 matching points).

We can see that the use of disparity maps results in a reduction of the error when illumination is different from that in which the training set was obtained. Additionally, the use of disparity result in a reduction of the variance of the error. This means that the error when using disparity is smaller in (almost) all test points (and not only in few of them).

The computation of the features for a given disparity image takes less than 0.1 seconds (also in a Pentium 4 at 2.0 GHz). Therefore, the use of disparity maps for appearance-localization is completely feasible.

VIII. CONCLUSIONS

In this paper, we have introduced the use of sparse disparity maps to increase the robustness of appearancebased localization to changes in illumination.

To compress sparse disparity maps to a reduced set of features, we have to deal with the problem of missing values. We have presented a novel EM-based algorithm to extract the principal components of a set of data that is more efficient in the way in which it deals with missing values than previously existing methods.

After the dimensionality reduction, it remains the open problem of how to merge the information provided by disparity maps with that provided by intensity images. We have proposed a simple linear combination of the sensors models build separately for each type of sensor. In this way, noise in one sensors do not corrupt the information provided by the other sensor.

We have shown that, using features from disparity maps in addition to those obtained from intensity images, we can improve the quality of the sensor model when illumination conditions are different from those in which the training set is obtained. Thus, disparity maps are a good option to increase the robustness of appearancebased robot localization.

In the EM algorithm we have introduced, the value for missing pixels for a given disparity map are inferred using information from other disparity maps. In our future work, we would like to explore the possibility doing this using also information provided by the intensity images (close points with similar intensity value are likely to have similar disparity).

The main assumption behind our approach is the existence of a training set obtained off-line and densely sampled over the space where the robot is expected to move. In general, obtaining this training set is not a problem, but it would be desirable the robot to build it on-line. With this extension, our system would become a Concurrent Localization and Mapping (CLM) system. To achieve this improvement, we have to explore the use of incremental techniques for compressing the images that are obtained on-line.

Fig. 2 .

 2 Fig. 2. Plain image (top) and the corresponding disparity map (bottom). In the disparity map, light gray areas are missing values.

Fig. 3 .

 3 Fig. 3. Our robot with the stereo head used to take images and disparity maps.

Fig. 4 .

 4 Fig. 4. Error in positioning in three different illumination conditions, using only intensity images (dashed line) and intensity and disparity images (dotted line).

IX. ACKNOWLEDGMENTS

This work has been partially supported by the European (ITEA) project "Ambience: Context Aware Environments for Ambient Services".

We would like to thank Bas Terwijn for helping us to obtain the data necessary to perform the experiments reported in this paper.

X. REFERENCES


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

N Vlassis
email: vlassis@science.uva.nl

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

A Variational EM Algorithm for Large-Scale Mixture Modeling

Keywords: Gaussian mixture, EM algorithm, variational approximation, clustering, very large database

Mixture densities constitute a rich family of models that can be used in several data mining and machine learning applications, for instance, clustering. Although practical algorithms exist for learning such models from data, these algorithms typically do not scale very well with large datasets. Our approach, which builds on previous work by other authors, offers an acceleration of the EM algorithm for Gaussian mixtures by precomputing and storing sufficient statistics of the data in the nodes of a kd-tree. Contrary to other works, we obtain algorithms that strictly increase a lower bound on the data log-likelihood in every learning step. Experimental results illustrate the validity of our approach.

Introduction

Mixture modeling and clustering are common tasks in data analysis and occur in applications across many fields. Clustering is employed to find prototypical data items in a large data base or to determine whether new data matches previously seen data. Clustering and mixture modeling are also used in supervised classifier systems, where new data is compared with each cluster and each cluster defines a mapping to the class labels. Mixture modeling is used for density estimation in general with many different applications.

The Expectation-Maximization (EM) algorithm and the k-means algorithm are among the most popular algorithms for data clustering and mixture modeling [START_REF] Mclachlan | Finite Mixture Models[END_REF]. This is mainly due to the ease of implementation of these algorithms and the fact that they are guaranteed to improve the associated objective functions in every step of the algorithms without the need to set any parameters. However, since for both algorithms the run-time per iteration is linear in both the number of data items n and the number of clusters k, their applicability is limited in large scale applications with many data and many clusters.

Recently several authors [START_REF] Moore | Very fast EM-based mixture model clustering using multiresolution kd-trees[END_REF][START_REF] Moore | Accelerating exact kmeans algorithms with geometric reasoning[END_REF][START_REF] Kanungo | An efficient k-means clustering algorithm: Analysis and implementation[END_REF][START_REF] Alsabti | An efficient k-means clustering algorithm[END_REF] proposed speedups of these algorithms based on analyzing large chunks of data at once to save distance computations. The idea is to use geometrical reasoning to determine that for chunks of data a particular prototype is the closest (k-means) or the posterior on mixture components hardly varies in the chunk of data (EM for Mixture of Gaussians). Cached sufficient statistics can then be used to perform the update step of the algorithms.

In this paper we employ a variational EM approach [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] and show how such speedups for the EM algorithm can be guaranteed to increase a lower bound on the data log-likelihood in each step. We also derive closed form and efficiently computable optimal assignments of chunks of data to mixture components (E-step). Furthermore, the variational EM framework allows for arbitrary partitionings, where existing techniques were forced to use relatively fine partitionings of the data. Thus the proposed framework allows more freedom in designing new speedup algorithms.

The rest of this paper is organized as follows: in the next section we discuss Gaussian mixture modeling and in Section 3 we derive the corresponding EM algorithm. Section 4 shows how we can use locally shared responsibilities to speed-up the EM algorithm. Then, in Section 5 we compare our work with related work. In Section 6 we describe our experiments on speeding up Gaussian mixture learning and end with conclusions in Section 7.

Clustering as Gaussian mixture learning

The problem of data clustering can be viewed as a density estimation problem if we treat the data points as samples from a mixture density of k components (clusters). A typical case is to assume Gaussian components.

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 loglikelihood. 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) decompo-sitions

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 the 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 loglikelihood. 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)

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 Chunky EM algorithm is to assign equal responsibilities to chunks of data points that are nearby in the input space.

Consider a partitioning A 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 ∈ A 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∈A 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 = xi∈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| xi∈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 log-likelihood 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 [START_REF] Sproull | Refinements to nearest-neighbor searching in k-dimensional trees[END_REF] can be written 1 :

log p(x|s) A = 1 |A| xi∈A log p(x i |s) = - 1 2 [log |C s | + 1 |A| xi∈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) 1 We ignore the additive constantd 2 log(2π) which translates into a multiplicative constant in [START_REF] Sproull | Refinements to nearest-neighbor searching in k-dimensional trees[END_REF]. w.r.t. θ to zero we find the update equations:

p(s) = 1 n A∈A |A|q A (s) (12)
m s = 1 np(s) A∈A |A|q A (s) x A (13)
C s = 1 np(s) A∈A |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.

Using a kd-tree for partitioning

The algorithm presented in the previous section works for any partitioning, as long as sufficient statistics of the data have been stored in the corresponding cells. We note that the amount of computation needed to perform one iteration of the Chunky EM increases linearly with the number of cells. Moreover, it is easy to see that if we refine a specific partitioning and recompute the responsibilities q A (s) then F cannot decrease. In the limit holds F = L and the approximation bound is tight. Clearly, various trade-offs between the amount of computation and the tightness of the bound can be made, depending on the particular application.

A convenient structure for storing statistics in a way that permits the use of different partitionings in the course of the algorithm is a kd-tree [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 that defines a hierarchical decomposition of the data space into cells, each cell corresponding to a node of the tree. Each node is split into two children nodes by a hyperplane that cuts through the data points in the node. Typically, a node is represented by a hyper-rectangle and an axis-aligned hyperplane is used for splitting nodes. In our case we use a kdtree where each node is split along the bisector of the first principal component of the data in the node [START_REF] Sproull | Refinements to nearest-neighbor searching in k-dimensional trees[END_REF], leading to irregularly shaped cells. 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 these cached statistics can be efficiently done bottom-up by noticing that the statistics of a parent node are the sum of the statistics of its children nodes. Building the kd-tree and storing statistics in its nodes has cost O(n log n).

A particular expansion of the kd-tree corresponds to a specific partitioning A of the data space, where each cell A ∈ A corresponds to an outer node of the expanded tree. Further expanding the tree means refining the current partitioning, and in our implementations as heuristic to guide the tree expansion we employ a best-first 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 partitioning (number of outer nodes at any step), or sampled approximations of the difference between log-likelihood and F. Another possibility, explained in the next section, is to control the tree expansion based on the performance of the algorithm, i.e., whether a new partitioning improves the value of F from the previous partitioning.

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 kd-tree 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 ). In comparison, we use q A proportional to the average joint log-likelihood, c.f. [START_REF] Sproull | Refinements to nearest-neighbor searching in k-dimensional trees[END_REF]. If there is little variation in the posteriors within a node (and this is achieved by having very fine-grained partitionings), the approximation will hardly affect the update in the M-step, and therefore the M-step will probably still 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 involve 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 we provably increase in each EM 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 run EM even with very coarse-grained partitionings. As mentioned above and as demonstrated in the experiments below, by gradually refining the partitioning while running EM we can get close to the optima of the log-likelihood in relatively few EM steps.

Experimental verification

In this section we describe our experiments and present and discuss the obtained results.

Implementation

In the experiment we apllied our variational EM algorithm to learn a k-component mixture of Gaussians for the data. First, a kd-tree is constructed as in [START_REF] Sproull | Refinements to nearest-neighbor searching in k-dimensional trees[END_REF] and we cache in its nodes the data statistics bottomup, as explained above. The EM algorithm is started with an initial expansion of the tree to depth two. We keep this partitioning fixed and run the variational EM till convergence. Convergence is measured in terms of relative increase in F. We then refine the partitioning by expanding the node of the tree that leads to maximal increase in F. Then we run again the variational EM till convergence, refine the partitioning by expanding best-first a single node of the tree, etc. We stop the algorithm if F is hardly improved between two successive partitionings.

Note that a refining a particular partitioning and computing the new (shared and optimal) responsibilities, can be viewed as applying a single E step which justifies the use of the relative improvement of F as a convergence measure.

In all experiments we used artificially generated data sampled from a randomly initialized kcomponent Gaussian mixture in d dimensions with a component separation of c. A separation of c means [START_REF] Dasgupta | Learning mixtures of Gaussians[END_REF]:

∀ i =j : ||m i -m j || 2 ≥ c 2 •max {i,j} {trace(C i ), trace(C j )}. (15)

Difference with suboptimal shared responsibilities

As indicated by theory, the responsibilities computed from (10) are optimal w.r.t. maximizing F. Since θ is unaffected by the choice of responsibilities, we see from (4) that the sum of Kullback-Leibler divergences should be smaller when using (10) than when for example using p(s| x A ), as in [START_REF] Moore | Very fast EM-based mixture model clustering using multiresolution kd-trees[END_REF]. In this experiment, we investigate the difference in summed Kullback-Leibler divergences (or equivalently the difference in F) between using [START_REF] Sproull | Refinements to nearest-neighbor searching in k-dimensional trees[END_REF] and the method of [START_REF] Moore | Very fast EM-based mixture model clustering using multiresolution kd-trees[END_REF].

We generated 20 data sets of 5000 points from a random k = 20 component Gaussian mixture in d = 10 dimensions with a separation of c = 2. We then initialized a mixture using k-means and approximated the responsibilities for this mixture. Figure 1 shows the differences in the Kullback-Leibler divergence summed over all data points, for different partitioning sizes, averaged over the 20 data sets. Each unit on the horizontal axis corresponds to expanding the tree one level down, and hence doubling the number of nodes used in the partition.

The result confirms that our method indeed gives a better approximation, as predicted from theory. The difference is larger for rough partitionings, this is as expected since p(s| x A ) becomes a better approximation of the individual responsibilities as the partition gets finer.

Gaussian mixture learning: Chunky vs. regular EM

In the second experiment we compared our method to the regular EM algorithm in terms of speed and quality of the resulting mixture. We looked at how difference in performance and speed is influenced by the number of data points, dimensions, components, and the amount of separation.

The default data set consisted of 10,000 points drawn from a 10-component 3-separated Gaussian mixture in 2 dimensions. To compare log-likelihoods we also created a test set of 1000 points from the same mixture. We applied both algorithms to a mixture initialized using k-means. Figure 2 shows the (differences in) negative log-likelihood and speed between the algorithms and the generating mixture averaged over 20 data sets. Speed was measured using the total number of basic floating point operations needed for convergence.

The results show with respect to the run-times that (i) the speedup is at least linear in the number of data points (ii) the number of dimensions and components have a negative effect on the speedup and (iii) the amount of separation has a positive effect. In general the Chunky EM requires more EM iterations to converge but it is still faster than regular EM since the iterations themselves are executed much faster. The difference in the performance of the two algorithms w.r.t. log-likelihood becomes larger (our method performs worse) as the number of components and data dimension increase. However, the difference in loglikelihood between Chunky and regular EM is still small compared to the respective differences with the log-likelihood of the generating mixture, even using many components in high dimensional spaces.

Conclusions

We presented a variational EM algorithm that can be used to speed-up large-scale applications of EM to learn Mixtures of Gaussians. Our contributions over similar existing techniques are two-fold.

Firstly, we can show directly that the algorithm maximizes a lower bound on data log-likelihood and that the bound becomes tighter if we use finer partitionings. The algorithm finds mixture configurations near local optima of the log-likelihood surface which are comparable with those found by the regular EM algorithm, but considerably faster.

Secondly, because we have a convergent algorithm that maximizes a lower bound on data log-likelihood we can use any partitioning of the data. This allows us to use partitionings where the true posteriors differ heavily in each part, which might be needed when working on very large data sets and only limited computational resources are available. Another option is to start the algorithm with a rough partitioning (for which the EM iterations are performed very fast) and only refine the partitioning when needed because the algorithm converged with the rough partitioning.

Comparing our work with [START_REF] Moore | Very fast EM-based mixture model clustering using multiresolution kd-trees[END_REF], we conclude that with the same computational effort we can use the optimal shared responsibilities instead of the posterior at the node centroid. Furthermore, we obtained a provably convergent algorithm and have the freedom to use arbitrary partitionings of the data.

Current research involves integrating the Chunky EM with greedy Gaussian mixture learning [START_REF] Vlassis | A greedy EM algorithm for Gaussian mixture learning[END_REF][START_REF] Verbeek | Efficient greedy learning of Gaussian mixture models[END_REF]. In greedy mixture learning, the mixture is build component-wise, starting from one component. After the EM for the k-component mixture converged, a new mixture component is added and EM cycles are performed on the resulting k + 1 component mixture, etc. Finally, we would like to remark that the proposed technique can directly be applied to other Gaussian mixture models like e.g. the Generative Topographic Mapping [START_REF] Bishop | GTM: The generative topographic mapping[END_REF] or the supervised mixture modeling [START_REF] Titsias | Shared kernel models for class conditional density estimation[END_REF].

Number of datapoints (x1000):

Figure 1 :

 1 Figure 1: Differences in Kullback-Leibler divergence as a function of partition size.

Figure 2 :

 2 Figure 2: Experimental results on Gaussian mixture learning. (Left) Speedup factor. (Right) Negative loglikelihoods at convergence: generating mixture (black), regular EM (light), Chunky EM (dark).

Acknowledgements

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


Laurent Hilico

Albane Douillet

Jean-Philippe Karr

Eric Tournié

A high transmission Faraday optical isolator in the 9.2 µm range

Keywords:

I. INTRODUCTION

Recent progress in quantum cascade laser (QCL) technology towards high power cw operation and off-the-shelf availability for a widening wavelength range has enlarged their application fields, for instance in molecular spectroscopy, lidar applications or THz imaging. Thanks to their wide tunability and integrability, QCLs are superseding alternative sources such as molecular gas lasers or lead-salt lasers [START_REF] Faist | Quantum cascade laser[END_REF][START_REF] Capasso | High-performance midinfrared quantum cascade lasers[END_REF][START_REF] Hugi | External cavity quantum cascade laser[END_REF]. One of the main drawbacks of QCLs is their extreme sensitivity to optical feedback that may prevent stable operation, hence the need for high optical isolation, especially for high-resolution spectroscopy applications in which precise frequency control is required.

QCLs are linearly polarized sources. A simple way to achieve optical isolation is to use a wire grid polarizer and a quarter-wave plate [START_REF] Mansfield | A lossless, passive isolator for optically pumped far-infrared lasers[END_REF]. However, the price to pay is rather high since the laser beam has to keep a circular polarization: any modification of the feedback beam polarization reduces the isolation ratio. The alternative solution relies on the Faraday effect, and was studied between the 60's and the 90's by several authors in the context of CO 2 laser applications. Mid-infrared Faraday isolation has been demonstrated using several Faraday media and different processes such as free carrier contribution at room temperature or interband and free carrier spin contributions at cryogenic temperatures [START_REF] Dennis | A 10.6-micron four-port circulator using free carrier rotation in InSb[END_REF][START_REF] Boord | Far-infrared radiation isolator[END_REF][START_REF] Jacobs | Faraday rotation optical isolator for 10.6-µm radiation[END_REF][START_REF] Tomasetta | 100 W average power 10.6 µm isolator based on the interband Faraday effect in InSb[END_REF][START_REF] Aggarwal | Faraday rotation in the 10 µm range in InSb at liquid-helium temperature[END_REF][START_REF] Carlisle | An optical isolator for mid-infrared diode lasers[END_REF][START_REF] Klein | Power handling capability of Faraday rotation isolators for CO2 laser radar[END_REF]. Isolation ratios up to 30 dB and insertion losses of 1.5 dB (71 % transmission) have been reported. n-doped InSb, which benefits from a strong free carrier Faraday effect, appears to be the most promising material in the 9-10 µm range [START_REF] Boord | Far-infrared radiation isolator[END_REF]. Nevertheless, it also suffers from rather large optical absorption and low heat conductivity, which make it unsuited for isolation of high-power CO 2 lasers. For this reason, development efforts were stopped, and currently no mid-IR Faraday isolator is commercially available.

Our initial motivation to fabricate such a device is the experiment developed in the Paris group, which we briefly describe here to illustrate an application requiring both a high isolation ratio and linearly polarized light. The experiment aims at a high-precision measurement of a Doppler-free two-photon vibrational transition frequency in the H + 2 molecular ion in the 9.2 µm range in order to obtain a new determination of the proton-to-electron mass ratio [START_REF] Roth | Precision spectroscopy of molecular hydrogen ions: towards frequency metrology of particle masses[END_REF][START_REF] Korobov | Relativistic corrections of mα 6 (m/M) order to the hyperfine structure of the H + 2 molecular ion[END_REF]. We developed a cw QCL source with a frequency control at the kHz level by a phase lock against a stabilized CO 2 laser [START_REF] Bielsa | Narrow-line phase-locked quantum cascade laser in the 9.2 µm range[END_REF][START_REF] Bielsa | HCOOH high-resolution spectroscopy in the 9.18 µm region[END_REF]. The QCL beam (≈54 mW) is mode-matched to a high finesse (≈1000) Fabry-Perot cavity, that is required both to ensure a Doppler-free geometry with perfectly counterpropagating beams, and to enhance the transition rate (two-photon ro-vibrational transitions in H + 2 being very weak). Optical feedback from the high-finesse cavity is a problem in such an experiment; in the initial setup, isolation was obtained by combining a quarter-wave plate and polarizer with an acousto-optic modulator providing a 6 dB additional isolation through its strongly polarization-dependent efficiency [START_REF] Karr | Vibrational spectroscopy of H + 2 : Hyperfine structure of two-photon transitions[END_REF]. The total isolation ratio, slightly over 30 dB, was just sufficient to achieve stable operation of the QCL. However, probing the transitions with a circularly polarized beam has two important drawbacks: (i) two-photon transition rates in H + 2 are almost an order of magnitude smaller with respect to the linear polarization case [START_REF] Karr | Vibrational spectroscopy of H + 2 : Hyperfine structure of two-photon transitions[END_REF], and (ii) transition frequencies are much more sensitive (by 4-5 orders of magnitude) to the Zeeman shift due to the |∆M | = 2 selection rule [START_REF] Karr | Vibrational spectroscopy of H + 2 : precise evaluation of the Zeeman effect[END_REF], requiring a high level of magnetic field control.

In this paper, we report the performances of the device we have fabricated, which is based on the same material (n-doped InSb) as most earlier realizations, but takes advantage of improvements in wafer quality and permanent magnets. After a brief review of Faraday effect in InSb, we describe the isolator setup and our experimental results. In the last section, we compare the performances of our device with published data and discuss the feasibility of a double-stage mid-infrared Faraday isolator.

II. FREE CARRIER FARADAY EFFECT IN INSB

At room temperature, infrared absorption and Faraday effect in n-doped InSb are dominated by the free electron contribution [START_REF] Boord | Far-infrared radiation isolator[END_REF] and can be described in the frame of the Drude model. The polarization rotation in a magnetic field B for a wafer thickness L is given by θ = V BL where

V = θ BL = µ 0 N q 3 λ 2 8π 2 n m * 2 c , (1)
is the Verdet constant [START_REF] Böer | A survey of semiconductor physics[END_REF], and the absorption coefficient α is

α = µ 0 N q 2 λ 2 4π 2 n c m * τ , (2)
where n = 4 is the material refractive index, m * is the carrier effective mass, N the carrier density and τ the effective carrier relaxation time. Equation (1) (resp. ( 2)

) are valid if the conditions ω r ≪ ω c ≪ ω ≤ ω p ≪ nω (resp. ω r ≪ ω)
are satisfied, where ω c is the cyclotron frequency, ω the optical field frequency, ω p the plasma frequency, and ω r = 1/τ . Note that both Faraday rotation and absorption follow a λ 2 law. The relevant figure of merit when assessing the performances of a material for optical isolation is the Faraday rotation F m per dB of attenuation and per Tesla. When assessing device performances (with a given magnetic field), it is best to use the Faraday rotation per dB attenuation F d . They are given by

F m = θ αLB = qcτ 2µ 0 m * (3)
F d = F m B. (4)
In order to optimize F d , one should have a B-field as strong as possible, and semiconductors of low effective mass and low resistivity, as discussed in [START_REF] Boord | Far-infrared radiation isolator[END_REF]. The figures of merit appear to be independent of the carrier density N , but a more detailed approach allows to study the influence of donor concentration and temperature [START_REF] Boord | Far-infrared radiation isolator[END_REF]; for InSb at 300 K, the optimal value is around N ≈ 2.10 17 cm -3 . We measured N = 2.35 10 17 to 2.5 10 17 cm -3 in our samples (see Table I) by Hall effect measurements. Let us check that the validity condition of equations (1-2) is met in our experimental conditions. The carrier effective mass in InSb is m * ≈ 0.014 m e [19]; taking as carrier density N = 2.35 10 17 cm -3 , the plasma frequency is ω p = (µ 0 N q 2 c 2 /m * ) 1/2 ≈ 2π×37 THz. With B ≈ 1 T the cyclotron frequency is ω c ≈ 2π × 2 THz. The effective carrier relaxation time τ is ≈ 1.4 10 -12 s, so ω r ≈ 2π × 0.1 THz and the laser frequency is ω ≈ 2π × 32 THz. One has indeed ω r ≪ ω c ≪ ω ≤ ω p ≪ nω.

III. EXPERIMENTAL SETUP

The experimental setup is depicted in figure III. The longitudinal magnetic field is produced by a standard permanent magnet. The solid line in figure III shows the B-field profile. The maximum field is 1.16 T. Four high-quality Te n-doped InSb wafers from Wafer Technology were cut into 4x4 mm square plates, antireflection coated and glued on a 0.5 mm pitch hollow copper screw to adjust the wafer position in the magnet. The wafers come from different slices of the same InSb ingot and exhibit slightly different characteristics that are reported in Table I. They are labeled 18, 43, 43B and 51 A.

All intensity measurements were done using two power-meters, one being used as a reference to cancel out possible laser power drifts. The wafer reflection coefficient was measured at small incidence to separate the incident and reflected beams. The residual reflectivity is ≤ 3%, instead of 36% without coating. The transmission coefficient T at λ = 9.166µ m (measured without polarizers) lies between 74 and 80% at room temperature for all samples. For Faraday rotation measurements, one wire grid polarizer was installed at each end of the Faraday rotator. The polarization rotation was monitored by rotating the second polarizer so as to minimize the transmitted intensity.

IV. RESULTS

First, we measured the Faraday rotation angles and absorption coefficients of three 520 µm-thick samples at three different wavelengths, using a CO 2 laser oscillating on the 9R42 (9.166 µm), 9P42 (9.753 µm) and 10P38 (10.787 µm) lines. Figure IV shows that both follow a λ 2 law, with respective slopes of 0.0065 cm -1 (µm) -2 and 0.76 deg(µm) -2 . Comparison with equations (1-2) yields m * =0.026 m e and τ =1. 4 10 -12 s. This effective mass is larger than the value of [19], but consistent with the value of 0.022 given in [START_REF] Zawadski | Electron transport phenomena in small-gap semiconductors[END_REF].

Since polarization rotations of about 60 deg were measured at 9.166 µm, the wafers were thinned to 390 µm to obtain rotation by about 45 deg for the maximum B-field. Figure III shows the polarization rotation as a function of the wafer position z in the magnet. As expected, it follows the B-field profile, and the value of the Verdet constant (see eq. ( 1)) can be extracted from this data. All characteristics are reported in Table I, as well as the corresponding figure of merit, showing that wafer 51A has the best performances.

The minimum transmitted intensity observed in polarization rotation measurements gives the Faraday isolation ratio. We measured between 27 dB and 30 dB, a value limited by the polarizer extinction ratio and by the power meter accuracy, but not by B-field or wafer dopant concentration inhomogeneities. The measured wire grid polarizer transmission is 90 %, so the overall transmission of the isolator is 0.8 × 0.9 2 = 65 % corresponding to 1.9 dB insertion loss.

Since Faraday rotation decreases and absorption increases with temperature, a simple way to improve the performances consists in lowering the wafer temperature [START_REF] Boord | Far-infrared radiation isolator[END_REF]. The isolator temperature was varied using Peltier coolers. Table II give the maximum Faraday rotation and wafer transmission at 14, 24 and 34 • C for the best wafer (51A). Wafer transmissions over 80% are achieved. The temperature dependance of the Faraday rotation angle is 1.

L N θmax V R T α F d wafer µm cm -3 deg. deg. T -1 mm -1 % % cm -1 deg.T -
5 deg K -1 .
This also shows that it is necessary to stabilize the isolator temperature to better than 1 K to ensure long-term stability of the isolation ratio in the 30 dB range. We checked that this optical isolator (operated at 16 • C to avoid water condensation) could be used for intracavity two-photon spectroscopy. The QCL beam was injected in the high finesse Fabry-Perot cavity (see Introduction), and locked on resonance. By carefully adjusting the polarizers, we obtained stable operation but only over a few minutes, indicating that the isolation ratio is only barely sufficient for this application.

As discussed in the Introduction, this result implies a significant increase of the transition probabilities in H + 2 twophoton spectroscopy. For precise comparison with the formerly used polarizer/quarter-wave plate system, insertion losses have to be taken into account. When the Faraday isolator is used, an additional half-wave plate of transmission 98 % is required to adjust the polarization, so that the overall transmission is T = 0.8 (wafer) ×0.9 2 (polarizers) ×0.98 = 0.635. The transmission of the former system is 0.9 (polarizer) ×0.98 (quarter-wave plate) = 0.88. Since the intensity of the considered two-photon lines is 8.5 larger in linear polarization [START_REF] Karr | Vibrational spectroscopy of H + 2 : Hyperfine structure of two-photon transitions[END_REF], and the transition rate being proportional to the square of the laser intensity, the net gain is a factor 8.5 × (0.635/0.88) 2 = 4.4.

V. DISCUSSION

Table III summarizes published characteristics of Faraday isolators in the 9-10 µm range and at 5.4 µm together with the present results. Note that the measurements reported in [START_REF] Dennis | A 10.6-micron four-port circulator using free carrier rotation in InSb[END_REF] were performed without antireflection coating on the wafer, leading to inaccurate absorption coefficient and figures of merit. The cryogenic temperature isolator of Ref. [START_REF] Tomasetta | 100 W average power 10.6 µm isolator based on the interband Faraday effect in InSb[END_REF] clearly exhibit the best figures of merit, but a rather low isolation ratio, with the drawback of a complicated setup. Due to the high B-field and good wafer quality, our isolator shows improved performance as compared to previously published ones at room temperature. The present results could be improved further in several ways. The simplest one is to use commercially available high-quality wire grid polarizers with 98% transmission and 1:500 extinction ratio. With such polarizers, a 27 dB isolation ratio is expected, and the insertion losses could be reduced to 1.1 dB.

Material N λ T 0 B isolation α Wafer ins. Fm F d (cm -3 ) (µm) (K) (T) (dB) (cm -
The Faraday isolator figure of merit F d can also be increased using higher B-fields and thinner wafers. Similar simple permanent magnet configurations giving 1.7 T and 2.1 T fields have been patented [21] or reported [START_REF] Trénec | Permanent magnets for Faraday rotators inspired by the design of the magic sphere[END_REF][START_REF] Mukhin | 2.1 Tesla permanent-magnet Faraday isolator for subkilowatt average power lasers[END_REF]. With 2.1 T, the wafer thickness can be reduced to 215 µm, resulting in a 94.7% wafer transmission (92% including the residual reflection losses after coating) for a single stage isolator, i.e. a 0.6 dB insertion loss assuming high quality (98% transmission) polarizers.

Very high isolation ratios are usually achieved using a double stage Faraday rotator, with two wafers and three polarizers. In that case, it is crucial to lower polarizer and wafer losses, since with presently reported values the overall transmission of such a device would be only 0.9 3 × 0.8 2 = 0.467 (for example, in H + 2 two-photon spectroscopy the gain in transition probability is reduced to a factor of 2.4). Using a 2.1 T field and high quality polarizers, one can expect a 79% isolator transmission (1 dB insertion loss) with an isolation ratio exceeding 50 dB for a double stage isolator.

VI. CONCLUSION

We have reported the analysis of Faraday rotation in n-doped InSb wafers and the operation of an optical Faraday isolator with a 30 dB isolation ratio and 1.9 dB insertion loss at 9.166 µm. We have also discussed possible improvements and shown the feasibility of a room temperature double stage isolator which would ensure >50 dB isolation ratio with 1 dB insertion loss. Such a device significantly extends the range of QCL applications , especially in high-resolution spectroscopy.

The present results can be extended to the transparency range of InSb. The short wavelength limit is λ ≈ 7.5 µm due to the 0.17 eV gap energy. The long wavelength limit is ≈ 30 µm due to the plasma frequency cut-off ω p /n. For shorter wavelengths in the 4-8 µm range, n-doped InAs could be used, but with higher absorption and reduced performances [START_REF] Boord | Far-infrared radiation isolator[END_REF].

FIG. 1 :FIG. 2 :

 12 FIG.1: Experimental setup. z is the wafer position with respect to the magnet center. The black boxes are power meters.

FIG. 3 :

 3 FIG.3: Maximum Faraday rotation (filled symbols) and absorption coefficient (open symbols) for three 520 µm-thick samples, versus the wavelength squared. The lines are adjustments by a proportional law with slopes of 0.76 deg(µm) -2 and 0.065 cm -1 (µm) -2 , respectively.

TABLE I :

 I Wafer characteristics: thickness L, carrier density N , polarization rotation maximum angle θmax, Verdet constant V , reflection and transmission coefficients R and T , absorption coefficient α and wafer figure of merit F d , for λ=9.166 µm.

				1
	43 390 2.5 10 17 44	97.3	3.1 78.4 4.63	210
	18 520 2.5 10 17 59	97.8	3.1 74.0 4.6	213
	43B 390 2.35 10 17 48	106	2.8 76.1 5.54	191
	51A 390 2.35 10 17 49	108	3.7 79.8 3.86	280

TABLE II :

 II Maximum polarization rotation and 51A wafer transmission versus the temperature T0.

	T0 ( • C) θmax T (%)
	14	50	80.4
	24	48.5 78.5
	34	47	76.8

TABLE III :

 III Comparison of published mid-infrared Faraday isolator performances.

Acknowledgments

The authors thank the training students A. Besse, A. Orgogozo and Ph. Genevaux for their participation in the experiments, as well as P. Grech for wafer polishing. This work was supported by the CNRS PEPS "métrologie du futur" action.