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

N Vlassis
email: vlassis@science.uva.nl

B Kröse
email: krose@science.uva.nl

Locally Linear Generative Topographic Mapping

We propose a method for non-linear data projection that combines Generative Topographic Mapping and Coordinated PCA. We extend the Generative Topographic Mapping by using more complex nodes in the network: each node provides a linear map between the data space and the latent space. The location of a node in the data space is given by a smooth nonlinear function of its location in the latent space.

Our model provides a piece-wise linear mapping between data and latent space, as opposed to the point-wise coupling of the Generative Topographic Mapping. We provide experimental results comparing this model with GTM.

Introduction

The Generative Topographic Mapping (GTM) (Bishop et al., 1998b) is a tool for non-linear data projection that has found many applications since its introduction. The idea is to fix a set of points in a latent space (which is used for visualization), these points are mapped through a Radial Basis Function Network [START_REF] Ripley | Pattern Recognition and Neural Networks[END_REF] into the data-space, where they are the means of the components of a Gaussian mixture. For projection, the data posterior probabilities on the mixture components are computed and the latent-space coordinates are averaged accordingly.

In this paper we extend the GTM by using more expressive covariance structure for the Gaussians in the mixture and, more importantly, by attaching to each mixture component a linear map between data and latent space.

The model may also be regarded as a restricted form of the Coordinated Principal Component Analysis (CPCA) model [START_REF] Verbeek | Coordinating Principal Component Analyzers[END_REF] (on its turn a restricted form of the coordinated factor analyzers model [START_REF] Roweis | Global coordination of local linear models[END_REF]). CPCA extends probabilistic Principal Component Analysis (PCA) [START_REF] Tipping | Mixtures of probabilistic principal component analysers[END_REF][START_REF] Roweis | EM Algorithms for PCA and SPCA[END_REF] by integrating the mixture of probabilistic PCA's into one nonlinear mapping between a data and a latent space. This is done by assigning a linear map from each PCA to a single, global, latent space. The objective function of CPCA sums data loglikelihood and a measure of how consistent the predictions of the latent coordinates are. With consistent we mean that PCA's that give high likelihood to a data point in the data space should give similar predictions on the (global) latent coordinate.

Here, we restrict the locations of the mixture components in the D-dimensional data space to be given by a smooth (non-linear) function of their location in the d-dimensional latent space (d D). The relation between the proposed model and the CPCA is analogue to the relation between the GTM and Kohonen's Self-Organizing Map [START_REF] Kohonen | Self-Organizing Maps[END_REF].

CPCA can be used for data visualization, however the model only aims at data loglikelihood maximization and at uni-modal distributions on the latent coordinate given a data point (consistency). What we would also like is uni-modal distributions in the other direction: the density on the data-space given a latent coordinate. By constraining the data-space centers of the Gaussian mixture components to be given by a smooth non-linear mapping of their latent space centers we safeguard CPCA against mixture components that are close in latent space but widely separated in the data space (which can cause multi-modal distributions from latent to data space).

The rest of the paper is organized as follows: in the following section we will discuss GTM and PCA and CPCA in some more detail.

Then, in Section 3 we present the new generative model and an EM-style learning algorithm is given in Section 4. Experimental results are presented and discussed in Section 5. We end the paper with a discussion and some conclusions in Section 6.

2 Generative Topographic Mapping, PCA and Coordinated PCA

Generative Topographic Mapping The Generative Topographic Mapping (GTM) (Bishop et al., 1998b) model is a tool for data visualization and data dimensionality reduction. Typically GTM is used for high dimensional numerical data. The data density is modeled by a mixture distribution that is parameterized in a way that enables low dimensional visualization of the data. The generative model of GTM is a mixture of k spherical Gaussians, all having the same variance σ 2 in all directions. The components of the Gaussian mixture are also called nodes, and play a similar role as the nodes in Kohonen's Self-Organizing Map [START_REF] Kohonen | Self-Organizing Maps[END_REF]. With each mixture component s, a fixed latent coordinate κ s is associated (typically the latent space is a two dimensional Euclidean space, and the κ s are chosen on a rectangular grid). In the latent space a set of m fixed (non-linear) smooth basis-functions φ i , . . . , φ m are defined, the vector φ(κ s ) denotes the response of the m basis functions on input κ s . The mean of mixture component s is given by Wφ(κ s ), where W is a D × m matrix. Hence, due to the smoothness of the basis-functions and because the multiplication with W is just a linear map, components with close latent coordinates will have close means in the data space.

Suitable parameters σ and W of GTM are typically found by employing a simple Expectation-Maximization (EM) algorithm. Like most EM algorithms, the algorithm can get stuck at locally optimal solutions, there is no guarantee to find globally optimal parameters

In order to visualize the high dimensional data in the latent space, for each data point the posterior distribution on the mixture components is computed. The latent coordinates of the nodes are then weighted according to the posterior distribution to give the latent coordinate for the data point.

Principal Component Analysis Principal Component Analysis (PCA) [START_REF] Jolliffe | Principal Component Analysis[END_REF]) is a widely known and used method for data visualization and dimensionality reduction. Assuming the data has zero mean, the data vectors {x n } are modeled as x n = Λg n + n . Here, x n is a D dimensional data vector, g n the corresponding latent coordinate, Λ a D×d dimensional matrix and n a residual vector. PCA uses the linear mapping Λ that minimizes sum of the squared norms of the residual vectors to map the data into a d dimensional representation.

The PCA linear map (matrix) Λ can also be found as the limiting case of σ → 0 of fitting a Gaussian density to the data where the covariance matrix is constrained to be of the form σ 2 I + ΛΛ The latter model is known as Probabilistic PCA [START_REF] Roweis | EM Algorithms for PCA and SPCA[END_REF].

Coordinated PCA Several Probabilistic PCA's can be combined in a mixture, by taking a weighted sum of the component densities [START_REF] Tipping | Mixtures of probabilistic principal component analysers[END_REF]. In the mixture case a set local PCA's is fitted to the data. Data can now be mapped in either one of the PCA models, however coordinates of the different PCA spaces cannot be related to each other. The coordinated PCA model [START_REF] Roweis | Global coordination of local linear models[END_REF][START_REF] Verbeek | Coordinating Principal Component Analyzers[END_REF] resolves this problem. Each PCA space is mapped linearly into a global low dimensional space. The global low dimensional space provides a single low dimensional coordinate system for the data. The question is which linear maps between the PCA's and the global latent space we would like.

It turns out that the data log-likelihood is invariant for how the different PCA spaces are coordinated, i.e. mapped into the global space. Therefore in [START_REF] Roweis | Global coordination of local linear models[END_REF][START_REF] Verbeek | Coordinating Principal Component Analyzers[END_REF] a penalty term is added to data log-likelihood to find a coordinated mixture of PCA's. For every data point x n we can compute the posterior probability on every PCA s. Furthermore, a prediction p(g n |s, x n ) on the global latent coordinate can be made using either PCA. The penalty term measures for every data point the ambiguity of the prediction on its global latent coordinate as given by p(g n |x n ) = s p(s|x n )p(g n |s, x n ). The measure of ambiguity is the smallest Kullback-Leibler divergence between p(g n |x n ) and any Gaussian density.

It turns out that adding this penalty term is exactly the same as using a variational EM algorithm [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]. Hence, it is relatively straightforward to derive EM algorithms to optimize the parameters of the model for this objective function.

In the next section we describe a generative model that combines the ideas of GTM and Coordinated PCA.

The Generative Model of Locally

Linear GTM

The generative model is best understood by starting in the latent space. First one of the k 'units' is chosen uniformly random, let's call it s. Then, a latent coordinate is drawn from p(g | s). Finally, an observation x is drawn from p(x | g, s). Using N (µ, Σ) to denote a Gaussian distribution with mean µ and covariance matrix Σ, the distributions, collectively parameterized by θ = {σ, ρ, α, W, Λ 1 , . . . , Λ k }, are:

p(s) = 1 k p(g | s) = N (κ s , α 2 ρσ 2 ) p(x | g, s) = N (Wφ(κ s ) + Λ s (g -κ s )/α, σ 2 I)
The marginal distribution on x given s reads:

p(x | s) = N (Wφ(κ s ), σ 2 (I + ρΛ s Λ s ))
Each mixture component p(x | s) is a Gaussian with variance σ 2 in all directions except those directions given by the columns of Λ s where the variance is ρ times larger. The form of the covariance matrix is that of probabilistic PCA, however here the vectors spanning the PCA subspace (given by the columns of Λ) all have the same norm.

The distribution on latent coordinates is a mixture of isotropic Gaussians (spherical covariance matrix) with equal mixing weights. The columns of the matrix Λ s span the PCA subspace of mixture component s, and provide the mapping for component s between the latent space and the data space. The parameter α is a magnification factor between latent and data space.

The generative model is almost the same as in [START_REF] Verbeek | Coordinating Principal Component Analyzers[END_REF], the only difference is that we constrain the data-space mean of the mixture component s to Wφ(κ s ) and keep the {κ s } fixed. The mapping is a linear sum, parameterized by the D × m matrix W, of m non-linear basis functions φ 1 , . . . , φ m of d inputs. The output of the function φ i (κ) is the i-th element of the vector φ(κ). In fact, we also add one constant basis function and d linear functions which just output one of the components of the input vector. These basis-functions can be used to capture linear trends in the data. In our model we consider the basis functions {φ i } to be fixed as well as the {κ s }, all other parameters are fitted to a given data set.

To compare, the generative model of GTM, which is a special case with ρ = 0 of the model used here, looks like:

p(s) = 1 k p(g|s) = δ(κ s ) p(x|s, g) = N (Wφ(κ s ), σ 2 I),
where δ(κ s ) denotes the distribution that puts all mass at κ s . This is similar to the model used in this paper, but (1) GTM uses a mixture of spherical Gaussians for the data space model and (2) GTM uses a mixture of delta peaks to model the density in the latent space.

Note that in our model, given a specific mixture component, the expected value of x depends linearly on g and vice versa, see Section 4.3. Whereas using GTM, given a mixture component s, there is no modeling of uncertainty in the latent coordinate g corresponding to a data point x, it is simply assumed that g = κ s .

Learning

The learning algorithm tries to achieve two goals. On the one hand we want a good fit on the observed data, i.e. high data log-likelihood. On the other hand, we want the mappings from data to latent space to give 'consistent' predictions, i.e. if a data point is well modeled by two mixture components in the data space, then these should give similar predictions on the corresponding latent coordinate of the data point.

In the following subsections we first discuss the objective function, then we describe how we initialize all the model parameters, and finally we discuss the learning algorithm.

The Objective Function

As shown by [START_REF] Roweis | Global coordination of local linear models[END_REF] the double objective can be pursued with a variational Expectation-Maximization (EM) algorithm. The free-energy interpretation of EM [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]) clarifies how we can use EM-like algorithms to optimize penalized loglikelihood objectives. The penalty term is the Kullback-Leibler (KL) divergence between the posterior on the hidden variables and a specific family of distributions, in our case uni-modal Gaussian distributions. The objective function sums data log-likelihood and the KL divergence which is a measure of how uni-modal the predictions on the latent coordinate are.

The objective Φ can be decomposed in two ways, given in equations ( 1) and ( 2) below, which can be associated respectively with the E and M step of EM. The first is the log-likelihood minus the KL-divergence, which is convenient for the E-step. The second decomposition isolates an entropy term independent of θ, which is convenient in the M-step.

Φ = - λ 2 W 2 + n log p(x n ; θ) - n D KL (Q n (g, s) p(g, s | x n ; θ)) (1) = - λ 2 W 2 + n E Qn log p(x n , g, s; θ) + n H(Q n (g, s)), (2)
where D KL and H are used respectively to denote Kullback-Leibler divergence and entropy [START_REF] Cover | Elements of Information Theory[END_REF]. The first term in both decompositions implements a Gaussian prior distribution over the elements of the matrix W, with inverse variance λ on each element. It makes sense to put a uniform prior on on the weights for the linear and constant basisfunctions. However, to keep notation simple, we assume equal prior on all weights here. The distributions Q n (g, s) over mixture components s and latent coordinates g are restricted to be independent over s and g, and of the form:

Q n (g, s) = q ns Q n (g) and Q n (g) = N (g n , Σ n ) It turns out that Φ is maximized w.r.t. Σ n if Σ n = v -1 I, where v = (ρ + 1)/(α 2 ρσ 2 )
is the inverse variance of p(g | x, s). Hence, in the following we assume this equality and do not use Σ n as a free parameter anymore. The {q ns } play the role of the 'responsibilities' in normal EM. Finally, g n is the expected latent coordinate for data point x n , if we approximate the true distribution on the latent coordinate with a uni-modal Gaussian. Hence, we might use the g n for visualization purposes.

Initialization

We initialize the model using Principal Component Analysis (PCA) [START_REF] Jolliffe | Principal Component Analysis[END_REF]. Suppose, without loss of generality, that the data has zero mean. We use a projection of the data to the two dimensional PCA subspace to initialize the locations g n of the distributions Q n (g). The columns of the loading matrices {Λ s } are initialized as the principal eigenvectors of the data covariance matrix. The variance σ 2 is initialized as the mean variance in directions outside the PCA subspace. The {κ s } are initialized on a square grid such that they have zero mean and the trace of their covariance matrix matches that of the covariance of the {g n }. The W is then taken as to map the {κ s } onto {Λ s κ s }. Then, we assign the data point n to the component with κ s closest to g n . Using the non-empty clusters, we compute the ρ and α. To initialize ρ, we use the mean value (over all non-empty clusters) of the total variance in the cluster in the data space minus Dσ 2 and divide the sum over dσ 2 . For α 2 we use the average of the mean variance in the cluster in the latent space and divide it over ρσ 2 . Finally, we can compute optimal q ns and start the learning algorithm with a M-step of the EM-algorithm described in the next section.

The inverse variance λ of the prior on the elements of W can be made larger to get smoother functions. In our experiments we used λ = 1.

Updating the Parameters

Using the notation:

E ns = 1 2σ 2 x ns x ns + v 2 g ns g ns -σ -2 α -1 x ns Λ s g ns +D log σ + d 2 log (ρ + 1),
x ns = x n -Wφ(κ s ), and g ns = g nκ s .

The objective (2) can be written as:

Φ = - λ 2 W 2 - ns q ns log q ns + E ns + c,
where c is some constant in the parameters. All the update equations below can be found by setting the derivative of the objective to zero for the different parameters. First, we give the Estep update equations:

q ns = e -Ens s e -E ns , g n = s q ns g n s ,
where for all N data points the expected latent coordinate given component s is:

g n s = E p(g|x,s) [g] = κ s + αρ ρ + 1 Λ s x ns .
Next, we give the M-step update equations. For Λ s we have to minimize:

n q ns x ns Λ s g ns .
This problem is known as the 'weighted Procrustes rotation' [START_REF] Cox | Multidimensional Scaling[END_REF]. Let

C = [ √ q 1s x 1s • • • √ q N s x N s ][ √ q 1s g 1s • • • √ q N s g N s ]
with SVD: C = ULΓ , (the g ns have been padded with zeros to form D-dimensional vectors). The optimal Λ s is given by the first d columns of UΓ . For W we find:

W = ns q ns (x n -α -1 s Λ s g ns )φ(κ s ) × ns q ns φ(κ s )φ(κ s ) + λσ 2 I -1
.

For the updates of the three variance parameters we use the following notation:

A = ns q ns x ns Λ s g ns , B = ns q ns x ns x ns , C = ns q ns g ns g ns .
Simultaneously setting the derivatives of Φ w.r.t. α, ρ and σ to zero gives:

ρ + 1 = (D -d)A 2 d(BC -A 2 ) , α = C(ρ + 1) Aρ , σ 2 = B -α -1 A N D

Experimental Illustration

In this section we provide experimental results comparing the presented model with GTM. The first two experiments use artificially generated data. In the last experiment we use chemical chromatogram data from the Unilever research department.

Artificial data The first data set consisted of points in IR 2 drawn along the sine function. We used GTM (20 units, 5 basis-functions) and our new model (5 units, 3 basis-functions) to map the data to a one dimensional latent space. The fitted models and the latent space representations are given in Figure 1. We see that for GTM the distribution of the latent coordinates clusters around the mixture component locations. For our model, while using 4 times fewer mixture components, the latent coordinates do not show this clustering due to the piece-wise linear maps.

Next, we generated a data set in IR 3 by first drawing uniformly random from the unit square in IR 2 , then to every data point x n we then added a third coordinate x

(3) n with:

x (3) n = (x (1) n -1/2) 2 + 2(x (2) n -1/2) 3 + ,
where is a random variable with distribution N ( ; 0, 100 -1 ). In Figure 2 we show the projections found by our method and GTM. Again we see that GTM concentrates the latent coordinates around the locations of the units, the projection 'clumps' around the units. Whereas our method gives a projection that matches the original distribution on the unit square better.

Chromatogram data

We demonstrate the performance of the proposed method and GTM on a problem from high-throughput screening of Maillard reactions, which involves the heating of sugars and amino acids, giving rise to complex mixtures of volatile flavor compounds.

The reaction product mixtures are characterized through fast gas chromatography. At Unilever R&D labs, it was found that a single GTM plot of the chromatographic data, based on two latent variables, is equally informative as the screening of many score plots from various PCA combinations. The data was preprocessed by mapping it with PCA onto a 20 dimensional linear subspace of the original 2800 dimensional space. We call an ingredient separable if all the measurements on compounds which contain, or do not contain, the ingredient form a cluster in the projected data. In our experiment we found that most ingredients which were separable using GTM were also separable with our new projection method. In the experiments GTM used 400 nodes and 16 basis-functions, our model used only 36 units and 16 basis functions.

Discussion and Conclusions

Discussion Constraining the means of the mixture components to the manifold spanned by a smooth mapping ensures that components that are close in latent space will also be close in the data space. This makes the constrained version of Coordinated PCA more suitable for visualization purposes.

Comparing the presented model with GTM, we observe that we extended the point-wise cou-pling between data and latent space to a locally linear coupling. This removes the effect of 'clumpy' projections, as discussed in the previous section.

Several variations on the presented model are possible. For example, we might model the variances per mixture component and make the {κ s } free parameters. However, here we wanted to keep the model close to GTM. As for the last option, the optimization w.r.t. the {κ s } is a somewhat more involved due to the nonlinear mappings φ. However, using local linear approximations of the non-linear functions we could do the optimization.

Note that the component covariance structure is not aligned with the non-linear function. However, in practice we do expect the local subspaces to be more or less aligned with the manifold due to the coordination (see Figure 1). To speed-up the algorithm, we might set the {Λ s } as manifold aligned and check whether this increases the objective function. We only need re-compute the {Λ s } with SVD if taking them manifold aligned does not increase the objective function.

In (Bishop et al., 1998a) an extension of GTM is discussed where the the covariance matrix of p(x|s) is taken to be of the form σ 2 I + ηDD , where the d columns of the D × d matrix D contain the partial derivatives of the RBF network. 1 This generative model has a similar form of the covariance matrix in the data space as we use, however in the latent space it uses the same discrete distribution as GTM. Therefore, the extension described in (Bishop et al., 1998a), just like normal GTM, doesn't have the linear dependence of the expected value of g on x given s.

Conclusions

We proposed a generative model and corresponding learning algorithm for visualization and projection of highdimensional data. The model integrates the GTM and Coordinated PCA.

Experimental results show that the model does not suffer from the 'clumpy' projections of GTM. Furthermore, in the latent-space we have a true density function where GTM only provides a discrete 'spike' distribution on the latent-space. Also, the proposed model provides 1 The M step of the proposed EM algorithm is not exact, due to the dependence of D on W. a density p(x | g), whereas for GTM this distribution is not clearly defined.

Finally, the proposed model provides a measure of uncertainty on the found latent coordinates summarized in the distributions Q n (g).

Figure 1 :

 1 Figure 1: Top left panel: GTM fit (disks give means of mixture components, dashed line the ordering of the components in latent space) on data (dots). Top right panel: our new model (the sticks depict vectors √ ρσΛ s used for the linear maps). Bottom panel: (top) latent coordinates for our method, and GTM (middle 20 nodes, bottom 10 nodes). Both sets of latent points have the same variance.

Figure 2 :

 2 Figure 2: Projections of (from left to right) the proposed model, GTM and original 2D coordinates (which are rotated).

Figure 3 :

 3 Figure 3: Projections of all the data for the proposed model (top) and GTM (bottom). From left to right we highlight compounds with one of three ingredient with boxes.

Acknowledgments

We would like to thank Dick de Ridder and Sam Roweis for fruitful discussions. We would also like to thank Claire Boucon and Sijmen de Jong (Unilever R&D Vlaardingen, The Netherlands) for providing the Maillard reaction data and experimental results. This research is supported by the Technology Foundation STW (project nr. AIF4997) applied science division of NWO and the technology program of the Dutch Ministry of Economic Affairs.


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

N Vlassis

B Kröse

Intelligent Autonomous Systems

Keywords: dimension reduction, feature extraction, principal manifold, selforganization

, this problem can be expressed as a penalized likelihood optimization problem. We show that a restricted form of the Mixtures of Probabilistic Principal Component Analyzers model allows for an efficient EM-style algorithm. The Procrustes Rotation, a technique to match point configurations, turns out to give the optimal orientation of the patches in the global space. We also show how we can initialize the mappings from the patches to the global coordinates by learning a non-penalized density model first. Some experimental results are provided to illustrate the method.

Section 1 Introduction 1 1 Introduction

With increasing sensor capabilities, powerful feature extraction methods are becoming increasingly important. Consider a robot sensing its environment with a camera yielding a stream of 100 × 100 pixel images, i.e. a stream of 10.000 dimensional vectors if we regard the image as a pixel intensity vector. The observations made by the robot often have a much lower intrinsic dimensionality. If we assume a fixed environment, and a robot that can rotate around its axis and drive through a room, the intrinsic dimensionality is only three. Linear feature extraction techniques are able to do a fair compression of the signal by mapping it to a much lower dimensional space. However, only in very few special cases the manifold on which the signal is generated is a linear subspace of the sensor space. This clearly limits the use of linear techniques and suggests to use non-linear feature extraction techniques.

Mixtures of Factor Analyzers (MFA) [START_REF] Ghahramani | [END_REF] and Mixtures of Probabilistic Principal Component Analyzers (MPPCA) [START_REF] Tipping | [END_REF] can be used to model such non-linear data manifolds. These methods (in fact MPPCA is a special case of MFA) provide a mapping back and forth between latent and data-space. However, these mappings only have a local applicability and are not related, i.e the coordinate systems of neighboring facor analyzers might be completely differently oriented. We cannot combine the different latent spaces, the coordinates of a point in one subspace give no infotmation on the coordinates in a neighboring subspace.

Recently, a model was proposed that integrates the local linear models into a global latent space, allowing for mapping back and forth between the latent and the data-space. The idea is that there is a linear map for each factor analyzer between the the data-space and the global latent space. This penalized log-likelihood model proposed in [START_REF] Roweis | [END_REF], solved by an EM-like procedure, is discussed in the next section.

Here, we show how we can remove the iterative procedure from the M-step by simplifying the density model. Furthermore, we show how we can use an 'uncoordinated' mixture model to initialize the mappings from data to latent space, providing an alternative to using an external (unsupervised) method to initialize the latent coordinates for all the data points.

In the next section, we describe the density model that is used. Then, in Section 3 we show how this density model removes one of the iterative procedures in the algorithm given in [START_REF] Roweis | [END_REF]. Section 4 discusses how an 'uncoordinated' mixture model can be used to initialize the mappings between latent and data-space. Some experimental results are given in 5. We end with a discussion and some conclusions in 6.

The density model

To model the data density in the high dimensional space we use mixtures of a restricted type of Gaussian densities. The mixture is formed as a weighted sum of its component densities, the weight of each component is called its 'mixing weight' or 'prior'. The covariance matrices of the Gaussians are constrained to be of the form:

C = σ 2 (I D + ρΛΛ ), Λ Λ = I d , ρ > 0 (1)
where D and d are respectively the dimension of the high-dimensional/data-space and the lowdimensional/latent space. We use I d to denote the d-dimensional identity matrix. The columns of Λ, in factor analysis known as the loading matrix, are D-dimensional vectors spanning the subspace. Directions within the subspace have variance σ 2 (1 + ρ), other directions have σ 2 variance. This is as the Mixture of Probabilistic Principal Component Analyzers (MPPCA) model, with the difference that here we do not only have isotropic noise outside the subspaces but also also isotropic variance inside the subspaces. We use this density model to allow for convenient solutions later. In Appendix A we derive a Generalized EM algorithm to find maximum likelihood solutions for this model. The same model can be rephrased using hidden variables z, which we use to denote 'internal' coordinates of the subspaces. We scale the coordinates z such that:

p(z | s) = N (0, I d ), (2)
where N (µ, Σ) denotes the Gaussian density, centered at µ and with covariance matrix Σ. The internal coordinates allow us to clearly express the link to the global latent space, for which we denote coordinates with g. All mixture components/subspaces have their own orthogonal mapping to the global space, parameterized by a translation κ and a matrix A, i.e. p(g | z, s) = δ(κ s + A s z), where δ(•) denotes the distribution with mass 1 at the argument. We use p s to denote the mixing weight of mixture component/subspace s and µ s to denote its mean or location. The generative model then looks like:

p(x) = s p s N (µ s , σ 2 s (I D + ρ s Λ s Λ s )) = s p s z dz p(x | z, s)p(z | s) (3) p(x | z, s) = N (µ s + √ ρ s σ s Λ s z, σ 2 s I D ) (4) p(g) = s p s N (κ s , A s A s ) = s p s N (κ s , α 2 s σ 2 s ρ s I d ). ( 5
)
We put an extra constraint on the projection matrices:

A s = α s σ s √ ρ s R s , R s R s = I d , α s > 0. (6)
The matrix R s implements only rotations plus reflections due to the orthonormality constraint.

The generative model reads: first 'nature' picks a subspace according to the prior distribution {p s }. A subspace generates internal coordinates z , which on the one hand give rise to hidden global coordinates g (via the translation κ s and the rotation plus scaling A s ) and on the other hand to observable data x (via the translation µ s and the orthogonal mapping √ ρ s σ s Λ s ) plus some noise with covariance σ 2 I D . Figure 1 illustrates the model. Note that the model assumes that locally there is a linear correspondence between the data space and the latent space. Furthermore, note that the densities p(g | x, s) and p(x | g, s) are Gaussian densities and hence that p(g | x) and p(x | g) are mixtures of Gaussians. In the next section we discuss how this density model differs from [START_REF] Roweis | [END_REF].

A simplified algorithm due to simplified density model

The goal is, given observable data {x n }, to find a good density model in the data-space and mappings {A s , κ s } that give rise to 'consistent' estimates for the hidden {g n }. With consistent we mean that if a point x in the data-space is well modeled by two subspaces, then the corresponding estimates for its latent coordinate g should be close to each other, i.e. the subspaces should 'agree' on the corresponding g.

Objective Function

To measure the level of agreement, one can consider for all data points how uni-modal the distribution p(g | x) is. This idea was also used in [START_REF] Vlassis | [END_REF], there the goal was to find orhtogonal projections for supervised data, such that the projections from high to low dimension preserve the manifold structure. In [START_REF] Roweis | [END_REF] it is shown how the double objective

' d d d d d d d d d s s E C S = k G S = 2 Z Z X µ, Λ, ρ, σ κ, A Z S = 1 Figure 1:
The generative model: x and g are independent given s, z.

of likelihood and uni-modality can be implemented as a penalized log-likelihood optimization problem.

Let

Q(g | x n ) denote a Gaussian approximation of p(g | x n ) = s p(g | x n , s)p ns = s p(g, s | x n ), with p ns = p(s | x n ). We define: Q(g, s | x n ) = Q(s | x n )Q(g | x n ), (7)
where

Q(g | x n ) = N (g n , Σ n ) and Q(s | x n ) = q ns .
As a measure of uni-modality we can use a sum of Kullback-Leibler divergences:

ns g Q(g, s | x n ) log Q(g, s | x n ) p(g, s | x n ) dg = (8) n D KL ({q ns } {p ns }) + s q ns D ns , (9)
where

D ns = D KL (Q(g | x n ) p(g | x n , s))
. The total objective function, combining loglikelihood and the penalty term, then becomes:

Φ = n log p(x n ) -D KL ({q ns } {p ns }) - s q ns D ns (10) = ns g dg Q(g, s | x n )[ -log Q(g, s | x n ) + log p(x n , g, s)] (11)
The form of (11) shows that we can view this as a variational approach to fit the mixture model where we consider both s and g as hidden variables. We approximate the true p(g, s | x) with the 'simple' Q(g, s | x). The terms involving the mixture model p, represent the expected likelihood.

The penalty term makes the distributions

Q(g, s | x n ) mimmic the distributions p(g, s | x n ).
Our density model differs with that of [START_REF] Roweis | [END_REF] in two aspects:

1. isotropic noise model outside the subspaces (as opposed to diagonal covariance matrix)

2. isotropic variance in subspace (as opposed to general Gaussian).

As a result, orthogonal mappings from the subspaces to the latent space can be absorbed in the matrices Λ, since for orthogonal matrices R we have ΛR (ΛR ) = ΛΛ , i.e. the density model does not change if we replace Λ with ΛR . Also, using our density model it turns out that to optimize Φ with respect to Σ n , it should be of the form1 Σ n = β -1 n I d . Therefore we work with β n from now on.

Using our density model we can rewrite the objective function (up to some constants) as:

Φ = ns q ns - d 2 log β n -log q ns - e ns 2σ 2 s - v s 2 dβ -1 n + g ns g ns ρ s + 1 (12) -D log σ s + d 2 log v s ρ s + 1 + log p s , (13)
where we used the following abbreviations:

g ns = g n -κ s , x ns = x n -µ s , e ns = x ns -α -1 s Λ s R s g ns 2 , v s = ρ s + 1 σ 2 s ρ s α 2 s . (14)

Optimization

Next, we give an EM-style algorithm to optimize Φ, it is a simplified version of the algorithm provided in [START_REF] Roweis | [END_REF]. The simplifications are: (i) the iterative process to solve for the Λ s , A s is no longer needed; an exact update is possible and (ii) the E-step no longer involves matrix inversions. The same manner of computation is used: in the E-step, we compute the uni-modal distributions Q(s, g | x n ), parameterized by β n , g n and q ns . Let g n s = E p(g|xn,s) [g] denote the expected value of g given x n and s. We used the following identities in the E-step update equations:

g n s = κ s + R s Λ s x ns α s ρ s /(ρ s + 1), ( 15
)
D ns = v s 2 [dβ -1 n + g n -g n s 2 ] + d 2 [log β n -log v s ]. ( 16
)
The approximating distributions Q can be found by iterating the fixed-point equations given below. Due to the form of the update for q ns it seems to make sense to initialize the q ns at p ns .

β n = s q ns v s , g n = β -1 n s q ns v s g n s , q ns = p ns exp -D ns s p ns exp -D ns . ( 17
)
In the M-step, we update the parameters of the mixture model: p s , µ s , κ s , Λ s R s , σ s , ρ s . Again we use some compacting notation:

C s = n q ns g ns 2 , E s = n q ns e ns , G s = d n q ns β -1 n (18)
Then, the update equations are:

p s = n q ns ns q ns , κ s = n q ns g n n q ns , µ s = n q ns x n n q ns , (19)
α s = C s + G s n q ns (g ns R s Λ s x ns ) , ρ s = D(C s + G s ) d(α 2 s E s + G s ) , σ 2 s = E s + ρ -1 s α -2 s [C s + (ρ s + 1)G s ] (D + d) n q ns (20)
Section 4 Initialization from a data-space mixture 5 Note that the above equations require E s which in turn requires the Λ s R s via equation ( 14). To find R s Λ s we have to minimize:

n q ns e ns = n q ns x ns -α -1 s ΛR g ns 2 or equivalently: - n q ns g ns (R s Λ s )x ns . (21)
This problem is known as the 'weighted Procrustes rotation' [Cox and Cox, 1994]. Let

C = [ √ q 1s x 1s • • • √ q ns x ns ][ √ q 1s g 1s • • • √ q ns g ns ] ,
with SVD:

C = ULΓ , (22)
where the g ns have been padded with zeros to form D-dimensional vectors, then the optimal Λ s R s is given by the first d columns of U Γ .

4 Initialization from a data-space mixture

To initialize the optimization procedure described in the previous section, we can go two ways.

The two ways correspond to starting with an E-step or with an M-step. In [START_REF] Roweis | [END_REF] it is proposed to start with an M-step, requiring initial estimates for the global coordinates g n . The initial global coordinates are typically provided by an 'external' unsupervised procedures, examples are LLE and Isomap. However, such methods might suffer from bad time complexity scaling with increasing sample-size2 . In this section we introduce an alternative initialization, that starts the update procedure with an E-step. This requires we find a data-space mixture model first, we can use the GEM algorithm described in Appendix A. Note that the greedy initialization method described in [Verbeek et al., 2001a], that builds the mixture model component-wise in order to avoid problems due to 'unlucky' random initialization of the mixture, can be directly used for the mixtue model described here.

Fixing the data-space mixture

Next, we consider how we can simplify the objective (13), if we fix the data-space mixture model. In this case the problem reduces to find appropriate mappings {κ s , A s }. Since the data-space mixture is fixed the log-likelihood term is constant and the objective reduces to the penalty term (9). If the density model is fixed, we can reduce the number of free parameters in (9) by setting q ns = p ns . This simplifies (9) to:

ns p ns D ns . (23)
Furthermore, we replace the Gaussian Q(g | x n ) with a delta-peak distribution on its center g n , this removes the free parameters Σ n . This amounts to simplifying (23) further to a weighted sum of log-likelihoods:

n,s p ns log p(g n | x n , s) (24)
The posterior distribution p(z | x, s) on internal coordinates z given a data point x is a Gaussian and hence also the posterior p(g | x, s) is a Gaussian. The covariance of the posterior p(z | x, s) is I d /(ρ s + 1) and to get the covariance for g we only have to multiply this with the square of the scaling factor of A s which yields:

I d α 2 s σ 2 s ρ s /(ρ s + 1) = I d v -1 s

Optimization and Initialization

Global maximization of ( 24) over all parameters is hard. Therefore, we maximize by alternating maximization over the global coordinate point estimates {g n }, and then over the patch locations and rotations {κ s , R s }. This leads to a local maximum of the objective. Skipping some constants, we rewrite (24) as:

- 1 2 n,s p ns [2d log α + v s g n s -g n 2 ]. (25)
In Appendix B we show how to find locally optimal parameters. Note that if the mixture in the data-space is kept fixed, then all computations involve d dimensional vectors and matrices, as expected.

Finally, to initialize the mappings we take component-wise approach, using the already initialized components to initialize the next. Note that due to the invariance w.r.t. global translations and rotations of the objective (24), we can simply initialize the first component with R = I d , α = 1, κ = 0. To initialize a new component s, we compute R s , α s , κ s which maximize:

n i∈F p ni p ns log p(g n | x n , s), ( 26
)
where F is the set of components for which we already fixed the mapping, i.e. the weights p ns are rescaled according to the already initialized components. Again, the solution is found in Appendix B. Note that these weights emphasize the data 'on the border' between the subspace s and the already initialized subspaces. In order to select the next subspace to initialize, we could for example pick s = arg max

s n p ns i∈F p ni /p s . (27)
This method somewhat resembles the method described in [Verbeek et al., 2001b]. There, also a similar density model is learned and then fixed to compute how the different local models can be combined into a global structure. The method presented here is more robust and applicable to any latent dimensionality.

Experimental results

In this section we discuss three experiments to illustrate the method. The first concerns artificially generated data. The other two experiments deal with mapping sets of images to low dimensional coordinates.

Artificial data-set

As a first experiment we used an artificially generated data-set, depicted in Figure 2. We learned a mixture model with 20 components and initialized the mappings to the latent space from the data-space mixture. For all data points we also have the coordinates on the surface, so we can inspect correlation coefficients with the discovered latent coordinates. The correlation coefficient between the first latent dimension and the surface coordinates are: 0.0494 and 0.9997. For the second latent dimension these are: 0.9961 and -0.0030.

Images of faces

For this experiment we captured gray valued images of a face with a webcam. Each image of 40 × 40 pixels is treated as a 1600-dimensional vector, with each coordinate describing the gray value of a specific pixel. Looking from left to rigth In the first experiment we used a set of 100 images of a face where there was only one degree of freedom in the data: the face was rotated from left to right when capturing the images. First we mapped the data to a 5 dimensional sub-space with PCA, capturing over 68% of the total variance in the data. Then we learned a coordinated mixture of principal component analyzers, the latent dimensionality was set to one and we used eight mixture components. To initialize, we used the method described in Section 4. The data-space mixture was learned with the greedy method of [Verbeek et al., 2001a], this method builds the mixture component per component to avoid bad initialization of the parameters. In figure 3 we show some of the used images, they are ordered according to the 1-d latent representation found by our method.

Two degrees of freedom Next, we consider a data set similar to the previous, except that the face now has a second degree of freedom, namely looking up and down. First we learned a coordinated mixture model with 1000 randomly selected images from a set of 2000 images of 40 × 40 pixels each. Again, we used a PCA projection, in this case to 22 dimensions, preserving over 70% of the variance in the data set. We used a latent dimensionality of two and 20 mixture components.

Here we initialized the coordinated mixture model by clamping the latent coordinates g n at coordinates found by Isomap and clamping the β n at small values for 50 iterations. The q ns were initialized at random values, and updated from the start. After the first 50 iterations, we also updated the β n and g n . The obtained coordinated mixture model was used to map the remaining 1000 images, the 'test-set' into the latent space. For each image x n we can compute p(g | x n ), as described in Section 3 this distribution can be approximated with a single Gaussian

Q n = arg min Q D KL (Q p(g | x n
)) with a certain mean and variance. We used the mean of Q n as the latent coordinate for an image.

In Figure 4 we show the latent coordinates for the test-set. We observe that a sensible 2dimensional parametrization is found in the training data set, that generalizes to the test set. To illustrate the discovered parametrization further, two examples of linear traversal of the latent space are given in Figure 5.

Robot navigation

Here we consider images obtained from a omni-directional camera3 mounted on a mobile robot. These images are warped to gray-scale panoramic images of 256 × 64 pixels, an example of such an image is given in Figure 6. The robot was placed at many locations in an office, while the orientation of the robot was kept fixed. At each location the position of the robot was recorded together with an image taken at that position.

First we consider data obtained from a rectangular area of an office, approximately 8 × 1 1 2 meters wide. The robot recorded images on a 10×10 centimeter grid, in the corridor we obtained 970 images. A uniformly random selected subset of 500 of those images were used to train the coordinated mixture model, with latent dimension two and 20 mixture components. Since we have 'supervised' data here, we can set and keep fixed the g n here at the recorded positions of the robot. In case of significant errors in the recorded positions we might only use them for initialization. Prior to training the model, we projected the images to a 15 dimensional space capturing over 70% of the variance in the total 16.384 dimensions. We mapped the test set to the latent space as described in the previous section. The covariance of the error between the true location and the found latent coordinates is almost diagonal with standard deviations of 13 and 17 centimeters in respectively the horizontal and vertical direction. The ellipse in Figure 7 illustrates the covariance structure, the axes of the ellipse are parallel to the eigenvectors of the covariance matrix, and their lengths equal to square root of the corresponding eigenvalues.

The same experiment was repeated with a data from a larger region, see the right pannel in Figure 7. Here the images were projected on a 50-dimensional space and 50 mixture components were used. From the total 2435 images 1700 were used for training and 735 for testing. In this experiment, the standard deviation of the error along the eigenvectors of the covariance matrix was 15 and 17 centimeters.

Although these results are good, considering that only linear models are used, we found it very hard to get good environment representation when not using the supervised position information.

6 Discussion and conclusions Discussion In the previous section we showed examples of application of the model to visual data. Both supervised and unsupervised data sets were used. The mains reason for applying this model to supervised data instead of other regression methods are the compact decription of the mapping between the high and low dimensional space, the probabilistic framework in which it is stated and the simplicity in terms of computational effort of the mapping.

One might consider applying this method to partially supervised data sets (only a limited number of measurements is provided together with a 'supervised' latent coordinate). This possibility forms a line of further research.

Another future line of research is the use of the decribed density model to get an idea of the An important issue, not addressed here, is that is many cases where we collect data from a systems with only few degrees of freedom we actually collect one or more sequences of data. If we assume that the system can vary its state only in continuous manner, these sequences should correspond to paths on the manifold of observable data. This fact might be exploited to find low dimensional embeddings of the manifold.

Conclusions

We showed how a special case of the density model used in [START_REF] Roweis | [END_REF] leads to a more efficient algorithm to coordinate probabilistic local linear descriptions of a data manifold. The M-step can be computed at once, the iterative procedure to find solutions for a Riccati equation is no longer needed. Furthermore, the update equations do not involve matrix inversions anymore. However, still d singular values of a D × D matrix have to be found.

We prosposed an alternative initialization scheme for the mapping from high to low dimension, that exploits the structure of a mixture model in the data space. We showed experimental results of this method discovering succesfully the structure in a data set of images of a face that rotates over its length axis. We showed how we can use the coordinated mixture model to map supervised data (robot navigation) and data for which we have uncertain latent coordinates (the 2 degrees of freedom facial images).

A GEM for constrained MPPCA

The constrained MPPCA model, characterized by equation [START_REF] Vlassis | Supervised dimension reduction of intrinsically low-dimensional data[END_REF], is parameterized by {p s , θ s }, the mixing weights p s and θ s = {µ s , σ s , ρ s , Λ s }. Let p ns denote the expectation that component s generated x n . The expected complete log-likelihood (expectation taken over the unobserved labels that indicate which component generated which data point) is then given by:

L = ns p(s | x n ) ln{π s p(x n ; θ i )}, (28)
where we used p to denote the density model with the new parameters. Optimization of the parameters is easy by using a Generalized EM (GEM) algorithm described below. First, we maximize (28) w.r. Next, we consider maximizing (28) w.r.t. the other parameters while keeping {p s } and { μs } fixed. This gives us parameter values that yield even higher values for (28) and hence we are guaranteed to obtain increased log-likelihood by the GEM algorithm. Before we turn to the derivation of the GEM algorithm we note that the likelihood under component s reads:

p(x n | s) = σ -D s (2π[ρ s + 1]) -d/2 exp -1 2σ 2 s [ x ns 2 - ρ s ρ s + 1 Λ x ns 2 ] ( 30
)
where the last equality follows from (33) and we ignored some additive constants in the equalities.

Setting the derivative of ( 41

where the λ si 's are as before. Note that σ2 s is intuitively interpreted as the mean variance outside the subspace. Similarly, (ρ s + 1) equals the mean variance inside the subspace divided over σ2

s .

B A variation on weighted Procrustes analysis

The quadratic form we need minimize with respect to {g n }, {κ s }, {R s }, {α s } is:

ns p ns [2d log α s + c s α -2 s κ s + α s R s z ns -g n 2 ], (44)
for constants c s , {z ns }, {p ns }. Solving for the g n by differentiation of (44) gives:

g n = s p ns c s α -2 s (κ s + α s R s z ns ) s p ns c s α -2 s . (45)
Then, keeping {g n } fixed we solve for the other parameters. For κ s , this gives: κ s = n p ns g n n p ns -n p ns α s R s z n n p ns = ḡsα s R s zs ,

i.e. the translation makes the weighted means equal, note that α s is as given by (49). To find the rotation R and scaling α, we assume that the weighted means are already equal, it is easy to show that this gives the same R and α. To realize this we set gns = g nḡs and zns = z nszs and use these to find the scaling and rotation. Then, for R s we need to maximize:

n p ns g ns R s zns , (47)
which is solved by the Procrustes rotation [Cox and Cox, 1994]. Let In order to prevent degenerate solutions which collapse all data into a single point, we need to put a constraint on the scale of the latent space. Either we fix {α s } or we could re-scale after each step such that n g n 2 = 1.

Figure 2 :

 2 Figure 2: An artificially generated data set, 1000 points in 3d on a 2d surface.

Figure 3 :

 3 Figure 3: Images selected at equally space intervals in the discovered latent space.

Figure 4 :

 4 Figure 4: The latent coordinates for the test set. Each circle represents a latent coordinate (its center) and the uncertainty β on it (the radius of the circle is β -1/2). For some coordinates we displayed the corresponding image at the coordinate.

Figure 5 :

 5 Figure 5: Examples of linear traversal of the latent space. The followed trajectories correspond to the arrows in Figure 4.

Figure 6 :

 6 Figure 6: An example of a panoramic image as used by the robot.

Figure 7 :

 7 Figure 7: Positions in the robot navigation task. The training locations are shown left, the found locations for the test set are shown in the middle with the original locations right. The ellipse (top right) shows the covariance structure in the error. The right pannel shows the locations and error covariance for the second experiment.

  is R s = ΓU . For α > 0, setting the derivative to zero gives:

  t. {p s } and { μs }, while keeping the other parameters fixed. If we write p(s | x n ) = p ns this gives:

	μs = n p ns x n n p ns	,	ps =	n p ns ns p ns	(29)

  ) w.r.t. ρs equal to zero, and doing the same for σs we obtain: ρs + 1 = n p ns x ns Λs Λ s x ns dσ 2 (D + dρ s ) = n p ns x ns x ns n p ns

	s	n p ns	,	σ2 s .	(42)
	Combining these we obtain:

Once we realize that the matrices Vs in[Roweis et al.,

2002] are of the form cI d with our density model, it can be seen easily by setting ∂Φ/∂Σn = 0 that Σn = β -1 I d .

Both LLE and Isomap scale in principle at least quadratic, due to the neighborhood graph construction.

In fact, a normal camera looks in vertical direction onto a parabolic mirror mounted above the camera.

Acknowledgements

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

A.1 Updating the loading matrices

The weighted covariance matrix for component s is defined as:

where we use x ns = x nμs . Let λ sj , u sj denote the eigenvalues respectively eigenvectors of S s sorted from large to small eigenvalues. We show that for σ s , ρ s > 0 setting Λs = [u s1 • • • u sd ] maximizes the expected log-likelihood. To maximize (28) with respect to Λs for given μs , consider the terms for Λs in (28):

where T = [VQ] is a D × D matrix with pairwise orthonormal columns and Q an arbitrary D × (Dd) matrix that fits the orthonormality constraints. The matrix I * denotes the matrix which is all zero except for the first d diagonal elements. From this is easy to see that the determinant of C is invariant for Λ:

So maximizing ( 32) is equivalent to minimizing: For σs , ρs > 0, this is equivalent to maximizing:

i.e. maximizing the weighted variance in the subspace spanned by the first d columns of T . This is solved by the first d eigenvectors of the weighted covariance matrix. Note that this solution for Λs is independent of the actual σs and ρs . So after computing { μs } and the {p s }, we use the first d eigenvectors of the weighted covariance matrix as Λs .

A.2 Variance inside and outside the subspace

Copyright IAS, 2000


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

Keywords: self-organizing map, Gaussian mixture, free-energy minimization, Expectation-Maximization algorithm, vector quantization

We present a variational Expectation-Maximization algorithm to learn probabilistic mixture models. The algorithm is similar to Kohonen's Self-Organizing Map algorithm and can be applied on any mixture model for which we can find a standard Expectation Maximization algorithm. We maximize the variational freeenergy which sums data log-likelihood and Kullback-Leibler divergence between the neighborhood function and the posterior distribution on the components, given data. We illustrate the algorithm with an application on word clustering.

Introduction

Kohonen's Self-Organizing Map (SOM) [START_REF] Kohonen | Self-Organizing Maps[END_REF] is a data analysis method that combines vector quantization with topology preservation. With each quantizer we associate a fixed location in a 'latent' space. The latent space is of much lower dimension (typically just two) than the data space. The SOM algorithm finds locations for the quantizers in the data space such that the summed squared distance from data to closest node is small and simultaneously topology is preserved. Topology preservation means that nearby components in latent space are also nearby in the data space. As a consequence, data points associated with nearby components in latent space are coming from similar locations in the data space. Hence, to check if data points are 'close' in the data space, we can check whether their associated components are close in the latent space. Topology preservation makes SOMs useful for data visualization and dimensionality reduction.

We present a constrained or variational Expectation-Maximization (EM) learning algorithm to learn probabilistic mixture models similar to SOM. The algorithm applies to a wide class of mixture models, in principle to any mixture model for whcih we can find a standard EM algorithm.

In the next section we discuss the Gaussian mixture model that will serve as the example mixture model throughout and also briefly discuss the EM algorithm. In section 3 we present our Self-Organizing Mixtures algorithm. We compare our algorithm to several closely related algorithms in Section 4. A brief example application is discussed in Section 5 and we present our conclusions in Section 6.

A simple generative model and EM

As generative model consider the Mixture of Gaussians (MoG), given by :

p(x) = 1 k k s=1 p(x | s),
for D-dimensional data, where p(x|s) is an isotropic Gaussian with inverse variance β and mean µ s . Given a data x = {x 1 , . . . , x N }, and an initial parameter vector θ = {β, µ 1 , . . . , µ k } the EM algorithm [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] finds a local maximizer θ of the log-likelihood L(θ). Assuming data are independent and identically distributed:

L(θ) = n log p(x n | θ) = n log 1 k s p(x n | s).
Learning is facilitated by introducing a set of N hidden variables. Each hidden variable indicates which of the k mixture components generated the corresponding data point. The EM algorithm maximizes the negative free-energy:

F (Q, θ) = E Q log p(x, s; θ) + H(Q) (1) = L(θ) -D KL (Q p(s | x; θ)), (2)
where we used H to denote the entropy of a distribution and D KL to denote the Kullback-Leibler (KL) divergence between two probability distributions. Q = n q n is a distribution over the hidden variables, q n gives a distribution on the generating mixture component s ∈ {1, . . . , k} for data point n. Note that, due to the non-negativity of the KL-divergence, for any Q, F (Q, θ) is a lower bound on L(θ).

The two forms in which we expanded F are associated with the M-step and the E-step of EM. In the M-step we change θ as to maximize, or at least increase, F (Q, θ). The first decomposition includes H(Q) which is a constant w.r.t. θ. In the E-step we maximize F w.r.t. Q, the second decomposition includes L(θ) which is constant w.r.t. Q. What remains is a KL divergence which is the effective objective function in the E-step of EM.

In standard applications of EM (e.g. mixture modeling) Q is unconstrained, which results in setting q n = p(s | x n ; θ) in the E-step since the non-negative KL divergence equals zero if and only if both arguments are equal. Therefore, after each E-step F (Q, θ) = L(θ). Variational methods are used when optimization over the unconstrained Q is intractable, Q is then restricted to a certain class Q of distributions allowing for tractable computations. Variational EM maximizes F instead of L, the objective sums log-likelihood and a penalty which is high if the true posterior is far from any member of Q.

Using free-energy for self-organization

By constraining choosing an appropriate class of distributions Q, we can also enforce a topological ordering between the mixture components, as we explain below. The approach is much like the one taken in [START_REF] Roweis | Global coordination of local linear models[END_REF][START_REF] Verbeek | Coordinating Principal Component Analyzers[END_REF], where constraints on Q are used to align the local coordinate systems of the components of a probabilistic mixture of factor analyzers.

We associate with each mixture component s a latent coordinate g s . It is convenient to take the components as located on a regular grid in the latent space. We set Q to be class of discretized isotropic Gaussians in the latent space, centered on either one of the k component locations g s and with a particular fixed variance. The distributions q n play a similar role as the neighborhood function in Kohonen's SOM [START_REF] Kohonen | Self-Organizing Maps[END_REF]. The q n are thus restricted to be contained in the finite set of distributions Q = {p 1 , . . . , p k }, where:

p r (s) = exp (-λ g s -g r 2 )
t exp (-λ g tg r 2 )

.

A small λ corresponds to a broad distribution (high entropy), and for large λ the distribution p r becomes more peaked (low entropy). By performing EM steps, we can never decrease the objective. The M-step is as usual and in the E-step we select for each data point x n the distribution q n ∈ Q that increases the objective the most. Note that in the E-step we can always stick to the q n of the previous EM round without changing (and hence definately not decreasing) the objective funciton. Since the objective function might have local optima and EM is only guaranteed to give locally optimal solutions, good initialization of the parameters of the mixture model is essential to finding a good solution. Analogous to the method of shrinking the extent of the neighborhood function with the SOM, we can start with a small λ (broad neighborhood function) and increase it iteratively until a desired value is reached. In implementations we started with λ such that the p r are close to uniform over the components, then we run the EM algorithm until convergence. Note that if the p r are almost uniform the initialization of θ becomes irrelevant. After convergence we set λ new ← ηλ old with η > 1 (typically η is close to unity). In order to initialize the EM procedure with λ new , we initialize θ with the value found in running EM with λ old .

Using q ns = q n (s), for our MoG case F can be rewritten as:

F (Q, θ) = 1 2 N D log β - ns q ns β x n -µ s 2 /2 + log q ns .
For fixed λ an EM algorithm can be derived by differentiation:

• E: Determine (by means of exhaustive or sparse search in Q, see below) for each x n the distribution p r * ∈ Q that maximizes F , set q n = p r * .

• M: Set: µ s = n q ns x n / n q ns and β = N D/ ns q ns x nµ s 2 .

The component r * on which p r * is centered for data point n, is referred to as the 'winner' for x n . The computational cost of the E-step is O(N k 2 ), a factor k slower than Kohonen's SOM and possibly prohibitively slow in large-scale applications. However, by restricting the search for a winner in the E-step to a limited number of candidate winners we can obtain an O(N k) algorithm. A straightforward choice is to use the l components with the largest joint likelihood p(x, s) as candidates, corresponding for our MoG to smallest Euclidean distance to the data point. If none of the candidates yields a higher value of F (Q, θ) we keep the winner of the previous step, in this way we are guaranteed never to decrease the objective in every step. We found l = 1 to work well and fast in practice, in this case we only check whether the winner from the previous round should be replaced with the closest node.

Let us consider why our algorithm yields topology preservation. The q n are by construction localized in the latent space. This implies that although the winner gets the most mass from q n , the neighbors of the winner also get some mass and far-away components get practically no mass. As a consequence nearby components in latent space are forced to model similar data.

Also, consider a topology preserving configuration of the mixture model. We can destroy the topology preservation if we permute the mixture components in the data space but keeping their locations in the latent space fixed. Note that L(θ) is the same for both configurations. Clearly, for the topology preserving configuration the mass of the posterior distribution is more or less localized in the latent space, where this will not be the case for the permuted model. Therefore, the KL divergence is expected to much smaller for the topology preserving configuration.

Discussion

Discussion. Our algorithm is very similar to Kohonen's SOM when applied on the example MoG. If we use only l = 1 candidate winner in the E-step the difference with Kohonen's winner selection is that we we only accept the closest node as a winner when it increases the energy and keep the previous winner otherwise. Our M-step coincides exactly with the update rule of the batch SOM: it puts component s at the weighted sum of the data, each data item weighted proportionally to the value of the neighborhood function (centered on the winner of the data point) at s. For an increasingly peaked neighborhood function both SOM and our algorithm aplied to the simple MoG reduce to the k-means vector quantization algorithm.

In [START_REF] Graepel | Self-organ. maps: generalizations and new optimiz. techniques[END_REF][START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF] another SOM-like algorithm is proposed with objective function:

- ns q ns [β r h sr x n -µ r /2 + log q ns ].
There, the neighborhood function, implemented by the h sr , is fixed, but the winner assignment is soft. Instead of selecting one 'winner' an unconstrained distribution over the components is used: the q ns . The β is used for annealing: for very small β the entropy term, with only one global optimum, becomes dominant, whereas for large β the quantization error, with many local optima, becomes dominant. By gradually increasing β more and more structure is added to the objective function.

Our work differs in several ways. We use one localized distribution in the latent space to constrain q n , as opposed to a mixture of localized distributions in the latent space. As a consequence the easy speed-up of our algorithm does not apply to the algorithms in [START_REF] Graepel | Self-organ. maps: generalizations and new optimiz. techniques[END_REF][START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF]. We use the neighborhood function width (controlled by λ) for annealing as opposed to using β. Note that although both λ and β can be used for annealing, only β can be optimized efficiently as a free parameter. Both our objective function and the one of [START_REF] Graepel | Self-organ. maps: generalizations and new optimiz. techniques[END_REF][START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF] can be interpreted as a log-likelihood plus a penalty term for non-topology preserving configurations. The interpretation provided in this work is simpler and applies trivially to any mixture model where for the interpretation of [START_REF] Graepel | Self-organ. maps: generalizations and new optimiz. techniques[END_REF][START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF] it is not clear whether it applies to any mixture model.

Another similar model is presented in [START_REF] Anouar | Probabilistic self organized map: application to classification[END_REF]. There, in the E-step we set:

q n = arg max r s p r (s)p(x n |s).
The M-step finds a new parameter vector that maximizes:

n log s q n (s)p(x n |s).
This algorithm does not optimize a single likelihood function, even if we keep the neighborhood function width fixed, since the mixing weights vary for each data item and change throughout the iterations of the algorithm. The algorithm has run-time O(nk 2 ), but can benefit from the same speed-up used here.

Another model, often presented as the probabilistic version of the self-organizing map, is the Generative Topographic Map (GTM) [START_REF] Bishop | GTM: The generative topographic mapping[END_REF]. However, the manner in which GTM achieves the topographic organization is quite different from those used in the SOM models. In GTM mixture components are parameterized by a linear combination of nonlinear functions of the locations of the components in the latent space. The parameters of the linear map are learned from data. The number of nonlinear basis functions and their smoothness are to be found by model complexity selection procedures.

Illustration

We illustrate our algortihm with an example where the mixture does not contain Gaussians as components but products of Bernoulli distributions, to underline that we are not limited to using Gaussian mixtures.

Figure 1 shows results of modeling word occurrence data on a part of the 20 newsgroup data set, word x n is shown at location g n = s p(s|x n )g s . Each of the 100 words is a data item with 16242 binary features indicating its occurrence in each of 16242 documents.

We learned a k = 25 component mixture model with our algorithm, where each mixture component is a product of Bernoulli distributions:

p(x n |s) = 16242 i=1 p(x (i) n |s) = 16242 i=1 p s,i x (i) n (1 -p s,i ) (1-x (i) n ) ,
where

x (i)
n denotes the i-th element of the vector x n . Due to the huge dimensionality of this data, the posteriors are rather peaked (low entropy), resulting in g n almost equal to one of the g s . The result is that the words locations g n are almost equal to one of the component locations g s . Smoothing the posteriors changes the plot from one in which the g n form 25 heaps of points on the g s , into one where the g n are drawn toward other similar heaps and the points are more spread. To 'smoothen' the representation we used p (s|x n ) ∝ p(s|x n ) α for α < 1 such that the entropies of the p (s|x n ) were close to two bits. 1The different occurance patterns of the words through the documents can be identified in the latent space. For example, the lower right corner is devoted to words that occur in computer related documents.

Conclusions

We presented a penalized log-likelihood probabilistic mixture modeling method, similar to Kohonen's SOM. The probabilistic formulation offers several benefits: (i) The method is directly applicable to any mixture model for i.i.d. data for which we can find a regular EM algorithm. (ii) The generative model allows seamless embedding in larger probabilistic systems. It is for example straightforward to learn mixtures of Self-Organizing maps. (iii) There is a direct and clear link to the mixture model data log-likelihood.

In current research we are investigating a similar method and in a sense opposite method. In the work presented here, we started with fixed locations in the latent space and fitted a mixture by maximizing the free energy. The opposite approach uses a fixed mixture and finds a corresponding latent representation that maximizes a similar free energy. The latter free-energy has single unique maximum as a function of the latent space representation which can be found by finding few eigenvectors of a matrix with edge size k.

Figure 1 :

 1 Figure 1: Self-organizing Bernoulli models, k = 25. Lower plot zooms dense area in upper plot.

It is not difficult to show that p ∝ p α is the best approximation in Kullback-Leibler sense to p such that H(p ) = c where α depends on c.

Acknowledgements

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

IAS reports

This report is in the series of IAS technical reports. The series editor is Stephan ten Hagen (stephanh@science.uva.nl). Within this series the following titles appeared: