Vasil Khalidov

Florence Forbes

Radu Horaud

Conjugate Mixture Models for Clustering Multimodal Data

Keywords: Mixture models, EM algorithm, multisensory fusion, audiovisual integration, Lipschitz continutity, global optimization Modèles de mélange, algorithm EM, fusion multisensorielle, intégration audiovisuelle, continuité Lipschitz, optimisation globale

The problem of multimodal clustering arises whenever the data are gathered with several physically different sensors. Observations from different modalities are not necessarily aligned in the sense that there is no obvious way to associate or to compare them in some common space. A solution may consist in considering multiple clustering tasks independently for each modality. The main difficulty with such an approach is to guarantee that the unimodal clusterings are mutually consistent. In this paper we show that multimodal clustering can be addressed within a novel framework, namely conjugate mixture models. These models exploit the explicit transformations that are often available between an unobserved parameter space (objects) and each one of the observation spaces (sensors). We formulate the problem as a likelihood maximization task and we derive the associated conjugate expectation-maximization algorithm. The convergence properties of the proposed algorithm are thouroughly investigated. Several local/global optimization techniques are proposed in order to increase its convergence speed. Two initialization strategies are proposed and compared. A consistent model-selection criterion is proposed. The algorithm and its variants are tested and evaluated within the task of 3D localization of several speakers using both auditory and visual data.

Introduction

The unsupervised clustering of multimodal data is a key capability whenever the goal is to group observations that are gathered using several physically different sensors. A typical example is the computational modeling of biological multisensory perception. This includes the issues of how a human detects objects that are both seen and touched [START_REF] Pouget | A computational perspective on the neural basis of multisensory spatial representations[END_REF][START_REF] Ernst | Humans integrate visual and haptic information in a statistically optimal fashion[END_REF], seen and heard [START_REF] Anastasio | Using Bayes' rule to model multisensory enhancement in the superior colliculus[END_REF][START_REF] King | The superior colliculus[END_REF][START_REF] King | Multisensory integration: strategies for synchronization[END_REF] or how a human localizes one source of sensory input in a natural environment in the presence of competing stimuli and of a variety of noise sources [START_REF] Haykin | The cocktail party problem[END_REF]. More generally, multisensory fusion (Hall and Mc-Mullen, 2004;[START_REF] Mitchell | Multi-Sensor Data Fusion[END_REF] is highly relevant in various other research domains, such as target tracking [START_REF] Smith | Approaches to multisensor data fusion in target tracking: A survey[END_REF] based on radar and sonar data [START_REF] Naus | Simultaneous localization of multiple emitters[END_REF][START_REF] Coiras | Rigid data association for shallow water surveys[END_REF], mobile robot localization with laser rangefinders and cameras [START_REF] Castellanos | Mobile Robot Localization and Map Building: A Multisensor Fusion Approach[END_REF], robot manipulation and object recognition using both tactile and visual data [START_REF] Allen | Integrating vision and touch for object recognition tasks[END_REF][START_REF] Joshi | Multisensor Fusion: A Minimal Representation Framework[END_REF], underwater navigation based on active sonar and underwater cameras [START_REF] Majumder | Multisensor data fusion for underwater navigation[END_REF], audio-visual speaker detection [START_REF] Beal | A graphical model for audiovisual object tracking[END_REF][START_REF] Perez | Data fusion for visual tracking with particles[END_REF][START_REF] Fisher | Speaker association with signal-level audiovisual fusion[END_REF], speech recognition [START_REF] Heckmann | Noise adaptive stream weighting in audio-visual speech recognition[END_REF][START_REF] Nefian | Dynamic Bayesian networks for audio-visual speech recognition[END_REF][START_REF] Shao | Stream weight estimation for multistream audio-visual speech recognition in a multispeaker environment[END_REF], and so forth.

When the data originates from a single object, finding the best estimates for the object's characteristics is usually referred to as a pure fusion task and it reduces to combining multisensor observations in some optimal way [START_REF] Beal | A graphical model for audiovisual object tracking[END_REF][START_REF] Kushal | Audio-visual speaker localization using graphical models[END_REF][START_REF] Smith | Approaches to multisensor data fusion in target tracking: A survey[END_REF]. For example, land and underwater robots fuse data from several sensors to build a 3D map of the ambient space irrespective of the number of objects present in the environment [START_REF] Castellanos | Mobile Robot Localization and Map Building: A Multisensor Fusion Approach[END_REF][START_REF] Majumder | Multisensor data fusion for underwater navigation[END_REF]. The problem is much more complex when several objects are present and when the task implies their detection, identification, and localization. In this case one has to consider two processes simultaneously: (i) segregation [START_REF] Fisher | Learning joint statistical models for audio-visual fusion segregation[END_REF] which assigns each observation either to an object or to an outlier category and (ii) estimation which computes the parameters of each object based on the group of observations that were assigned to that object. In other words, in addition to fusing observations from different sensors, multimodal analysis requires the assignment of each observation to one of the objects. This observation-to-object association problem can be cast into a probabilistic framework. Recent multisensor data fusion methods able to handle several objects are based on particle filters [START_REF] Checka | Multiple person and speaker activity tracking with a particle filter[END_REF][START_REF] Chen | Real-time speaker tracking using particle filter sensor fusion[END_REF][START_REF] Gatica-Perez | Audiovisual probabilistic tracking of multiple speakers in meetings[END_REF]. Notice, however, that the dimensionality of the parameter space grows exponentially with the number of objects, causing the number of required particles to increase dramatically and augmenting computational costs. A number of efficient sampling procedures were suggested [START_REF] Chen | Real-time speaker tracking using particle filter sensor fusion[END_REF][START_REF] Gatica-Perez | Audiovisual probabilistic tracking of multiple speakers in meetings[END_REF] to keep the problem tractable. Of course this is done at the cost of loss in model generality, and hence these attempts are strongly application-dependent. Another drawback of such models is that they cannot provide estimates of accuracy and importance of each modality with respect to each object. The sampling and distribution estimation are performed INRIA in the parameter space, but no statistics are gathered for the observation spaces. Recently [START_REF] Hospedales | Structure inference for Bayesian multisensory scene understanding[END_REF] extended the single-object model of [START_REF] Beal | A graphical model for audiovisual object tracking[END_REF] to multiple objects: several single-object models are incorporated into the multiple-object model and the number of objects is selected by an additional hidden node, which thus accounts for model selection. We remark that this method also suffers from exponential growth in the number of possible models.

In the case of unimodal data, the problems of grouping observations and of associating groups with objects can be cast into the framework of standard data clustering which can be solved using a variety of parametric or non-parametric techniques. The problem of clustering multimodal data raises the difficult question of how to group together observations that belong to different physical spaces with different dimensionalities, e.g., how to group visual data with auditory data? When the observations from two different modalities can be aligned pairwise, a natural solution is to consider the Cartesian product of two unimodal spaces. Unfortunately, such an alignment is not possible in most practical cases. Different sensors operate at different frequency rates and hence the number of observations gathered with one sensor can be quite different from the number of observations gathered with another sensor. Consequently, there is no obvious way to align the observations pairwise. Considering all possible pairs would result in a combinatorial blow-up and typically create abundance of erroneous observations corresponding to inconsistent solutions.

Alternatively, one may consider several unimodal clusterings, provided that the relationships between a common object space and several observation spaces can be explicitly specified. Multimodal clustering then results in a number of unimodal clusterings that are jointly governed by the same unknown parameters characterizing the object space.

The original contribution of this paper is to show how the problem of clustering multimodal data can be addressed within the framework of mixture models [START_REF] Mclachlan | Finite Mixture Models[END_REF]. We propose a variant of the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm (with discussion)[END_REF][START_REF] Mclachlan | The EM algorithm and extensions[END_REF] specifically designed to estimate object-space parameters that are indirectly observed in several sensor spaces. The convergence properties of the proposed algorithm are thoroughly investigated and several efficient implementations are described in detail. The proposed model is composed of a number of modality-specific mixtures. These mixtures are jointly governed by a set of common object-space parameters (which will be referred to as the tying parameters), thus insuring consistency between the sensory data and the object space being sensed. This is done using explicit transformations from the unobserved parameter space (object space) to each of the observed spaces (sensor spaces). Hence, the proposed model is able to deal with observations that live in spaces with different physical properties such as dimensionality, space metric, sensor sampling rate, etc. We believe that linking the object space with the sensor spaces based on object-space-to-sensor-space transformations has more discriminative power than existing multisensor fusion techniques and hence performs better in terms of multiple object identification and localization. To the best of our knowledge, there has been no attempt to use a generative model, such as ours, for the task of multimodal data interpretation.

In Section 2 we formally introduce the concept of conjugate mixture models. Standard Gaussian mixture models (GMM) are used to model the unimodal data. The parameters of these Gaussian mixtures are governed by the object parameters through a number of object-space-to-sensor-space transformations (one transformation for each sensing modality). Through the paper we will assume a very general class of transformations, namely non-linear Lipschitz continuous functions (see below). In Section 3 we cast the multimodal data clustering problem in the framework of maximum likelihood and we explicitly derive the expectation and maximization steps of the associated EM algorithm. While the E-step of the proposed algorithm is standard, the M-step implies non-linear optimization of the expected complete-data log-likelihood with respect to the object parameters. We investigate efficient local and global optimization methods. More specifically, in Section 4 we prove that, provided that the object-to-sensor functions as well as their first derivatives are Lipschitz continuous, the gradient of the expected complete-data log-likelihood is Lipschitz continuous as well. The immediate consequence is that a number of recently proposed optimization algorithms specifically designed to solve Lipschitzian global optimization problems can be used within the Mstep of the proposed algorithm [START_REF] Zhigljavsky | Stochastic Global Optimization[END_REF]. Several of these algorithms combine a local maximum search procedure with an initializing scheme to determine, at each iteration, good initial values from which the local search should be performed. This implies that the proposed EM algorithm has guaranteed convergence properties. Section 5 discusses several possible local search initialization schemes, leading to different convergence speeds. In Section 6 we propose and compare two possible strategies to initialize the EM algorithm. Section 7 is devoted to a consistent criterion to determine the number of objects. Section 8 illustrates the proposed method with the task of audiovisual object detection and localization using binocular vision and binaural hearing. Section 9 analyses in detail the performances of the proposed model under various practical conditions with both simulated and real data. Finally, Section 10 the paper and provides directions for future work.

Mixture Models for Multimodal Data

We consider N objects n = 1 . . . N . Each object n is characterized by a parameter vector of dimension d, denoted by s n ∈ S ⊆ R d . The set s = {s 1 , . . . , s n , . . . , s N } corresponds to the unknown tying parameters. The objects are observed with a number of physically different sensors. Although, for the sake of clarity, we will consider two modalities, generalization is straightforward. Therefore, the observed data consists of two sets of observations denoted respectively by f = {f 1 , . . . , f m , . . . , f M } and g = {g 1 , . . . , g k , . . . , g K } lying in two different observation spaces of dimensions r and p,

f m ∈ F ⊆ R r and g k ∈ G ⊆ R p .
One key ingredient of our approach is that we consider the transformations:

F : S → F G : S → G (1)
INRIA that map S respectively into the observation spaces F and G. These transformations are defined by the physical and geometric properties of the sensors and they are supposed to be known. We treat the general case when both F and G are non-linear.

An assignment variable is associated with each observation, thus indicating the object that generated the observation: A = {A 1 , . . . , A m , . . . , A M } and B = {B 1 , . . . , B k , . . . , B K }. Hence, the segregation process is cast into a hidden variable problem. The notation A m = n (resp. B k = n) means that the observation f m (resp. g k ) was generated by object n. In order to account for erroneous observations, an additional N + 1-th fictitious object is introduced to represent an outlier category. The notation A m = N + 1 (resp. B k = N + 1) means that f m (resp. g k ) is an outlier. Note that we will also use the following standard convention: upper case letters for random variables (A and B) and lower case letters for their realizations (a and b). The usual conditional independence assumption leads to:

P (f, g|a, b) = M m=1 P (f m |a m ) K k=1 P (g k |b k ).
(2)

In addition, all assignment variables are assumed to be independent, i.e.:

P (a, b) = M m=1 P (a m ) K k=1 P (b k ). (3)
As discussed in Section 10, more general cases could be considered. However, we focus on the independent case for it captures most of the features relevant to the conjugate clustering task and because more general dependence structures could be reduced to the independent case via the use of appropriate variational approximation techniques [START_REF] Jordan | An introduction to variational methods for graphical models[END_REF][START_REF] Celeux | EM procedures using mean field-like approximations for Markov model-based image segmentation[END_REF].

Next we define the following probability density functions, for all n = 1 . . . N, N + 1, for all f m ∈ F and for all g k ∈ G:

P F n (f m ) = P (f m |A m = n), (4)
P G n (g k ) = P (g k |B k = n). (5)
More specifically, the likelihoods for an observation to belong to an object n are Gaussian distributions whose means F(s n ) and G(s n ) correspond to the object's parameter vector s n mapped to the observations spaces by the transformations F and G:

P F n (f m ) = N (f m ; F(s n ), Σ n ), (6) P G n (g k ) = N (g k ; G(s n ), Γ n ), (7)
with:

N (f m ; F(s n ), Σ n ) = 1 (2π) r/2 |Σ n | 1/2 exp - 1 2 f m -F(s n ) 2 Σn , (8)
where the notation vw 2 Σ stands for the Mahalanobis distance (vw) ⊤ Σ -1 (vw) and ⊤ stands for the transpose of a matrix. The likelihoods of outliers are taken as two uniform distributions:

P F N +1 (f m ) = U(f m ; V ), (9)
P G N +1 (g k ) = U(g k ; U ), (10)
where V and U denote the respective support volumes. We also define the prior probabilities π = (π 1 , . . . , π n , . . . , π N +1 ) and λ = (λ 1 , . . . , λ n , . . . , λ N +1 ):

π n = P (A m = n), ∀m = 1 . . . M, (11)
λ n = P (B k = n), ∀k = 1 . . . K. ( 12
)
Therefore, f m and g k are distributed according to two (N + 1)-component mixture models, where each mixture is made of N Gaussian components and one uniform component:

P (f m ) = N n=1 π n N (f m ; F(s n ), Σ n ) + π N +1 U(f m ; V ), (13)
P (g k ) = N n=1 λ n N (g k ; G(s n ), Γ n ) + λ N +1 U(g k ; U ). (14)
The log-likelihood of the observed data can then be written as:

L(f, g, θ) = M m=1 log N n=1 π n N (f m ; F(s n ), Σ n ) + π N +1 U(f m ; V ) + + K k=1 log N n=1 λ n N (g k ; G(s n ), Γ n ) + λ N +1 U(g k ; U ) (15)
where:

θ = {π 1 , . . . , π N , π N +1 , λ 1 , . . . , λ N , λ N +1 , s 1 , . . . , s N , Σ 1 , . . . , Σ N , Γ 1 , . . . , Γ N } (16
) denotes the set of all unknown parameters to be estimated using a maximum likelihood principle.

The graphical representation of our conjugate mixture model is shown in Figure 1. We adopted the graphical notation introduced in [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF] to represent similar nodes in a more compact way: the M (resp. K) similar nodes are indicated with a plate. The two sensorial modalities are linked by the tying parameters s 1 , . . . s N shown in between the two plates.

INRIA M K λ 1 , . . . , λ N , λ N +1 Γ 1 , . . . , Γ N G k B k F m s 1 , . . . , s N Σ 1 , . . . , Σ N π 1 , . . . , π N , π N +1 A m

Generalized EM for Clustering Multimodal Data

Given the probabilistic model just described, we wish to determine the parameter vectors associated with the objects that generated observations in two different sensory spaces. It is well known that direct maximization of the observed-data loglikelihood ( 15) is difficult to achieve. The expectation-maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm (with discussion)[END_REF][START_REF] Mclachlan | The EM algorithm and extensions[END_REF] is a standard approach to maximize likelihood functions of type (15). It is based on the following representation, for two arbitrary values of the parameters θ and θ:

L(f, g, θ) = Q(θ, θ) + H(θ, θ), (17) with Q(θ, θ) = E[log P (f, g, A, B; θ) | f, g; θ], (18)
and

H(θ, θ) = -E[log P (A, B | f, g; θ)|f, g; θ], (19)
where the expectations are taken over the hidden variables A and B. Each iteration q of EM proceeds in two steps:

• Expectation. For the current values θ (q) of the parameters, compute the conditional expectation with respect to variables A and B:

Q(θ, θ (q) ) = a∈{1...N+1} M b∈{1...N+1} K P (a, b|f, g; θ (q) ) log P (f, g, a, b; θ) (20)
• Maximization. Update the parameter set θ (q) by maximizing (20) with respect to θ:

θ (q+1) = argmax θ Q(θ, θ (q) ) (21)
It is well known that the EM algorithm increases the target function L(f, g, θ) in (15), i.e., the sequence of estimates {θ (q) } q∈N satisfies L(f, g, θ (q+1) ) ≥ L(f, g, θ (q) ). Standard EM deals with the parameter estimation of a single mixture model, and a closed form solution for (21) exists in this case. When the maximization ( 21) is difficult to achieve, various generalizations of EM are proposed. The M step can be relaxed by requiring just an increase rather than an optimum. This yieds Generalized EM (GEM) procedures [START_REF] Mclachlan | The EM algorithm and extensions[END_REF] (see [START_REF] Boyles | On the convergence of EM algorithms[END_REF] for a result on the convergence of this class of algorithms). The GEM algorithm searches for some θ (q+1) such that Q(θ (q+1) , θ (q) ) ≥ Q(θ (q) , θ (q) ). Therefore it provides a sequence of estimates that still verifies the non-decreasing likelihood property although the convergence speed is likely to decrease. In the case of conjugate mixture models, we describe in more detail the specific forms of the E and M steps in the following sections.

The Expectation Step

Using (3)-( 12) the conditional expectation (20) can be decomposed as:

Q(θ, θ (q) ) = Q F (θ, θ (q) ) + Q G (θ, θ (q) ), (22)
with

Q F (θ, θ (q) ) = M m=1 N +1 n=1 α (q) mn log π n P (f m |A m = n; θ) , (23)
Q G (θ, θ (q) ) = K k=1 N +1 n=1 β (q) kn log λ n P (g k |B k = n; θ) , (24) where α (q)
mn and β

(q)
kn denote the posterior probabilities α q) ). Their expressions can be derived straightforwardly from Bayes' theorem, ∀n = 1 . . . N :

(q) mn = P (A m = n|f m ; θ (q) ) and β (q) kn = P (B k = n|g k ; θ (
α (q) mn = π (q) n N (f m ; F(s (q) n ), Σ (q) n ) N i=1 π (q) i N (f m ; F(s (q) i ), Σ (q) i ) + V -1 π (q) N +1 , (25)
β (q) kn = λ (q) n N (g k ; G(s (q) n ), Γ (q) n ) N i=1 λ (q) i N (g k ; G(s (q) i ), Γ (q) i ) + U -1 λ (q) N +1 . ( 26
) INRIA and α (q) m,N +1 = 1 - N n=1 α (q)
mn and β

(q) k,N +1 = 1 - N n=1 β (q)
kn . Using ( 6)-( 10) the expressions above further lead to:

Q F (θ, θ (q) ) = - 1 2 M m=1 N n=1 α (q) mn f m -F(s n ) 2 Σn + log((2π) r |Σ n |π -2 n ) - - 1 2 M m=1 α (q) m,N +1 log(V 2 π -2 N +1 ), (27)
Q G (θ, θ (q) ) = - 1 2 K k=1 N n=1 β (q) kn g k -G(s n ) 2 Γn + log((2π) p |Γ n |λ -2 n ) - - 1 2 K k=1 β (q) k,N +1 log(U 2 λ -2 N +1
).

(28)

The Maximization Step

In order to carry out the maximization (21) of the conditional expectation (20), its derivatives with respect to the model parameters are set to zero. This leads to the standard update expressions for priors, more specifically ∀n = 1, . . . , N + 1:

π (q+1) n = 1 M M m=1 α (q) mn , (29)
λ (q+1) n = 1 K K k=1 β (q)
kn .

(30)

The covariance matrices are governed by the tying parameters s (q+1) n ∈ S through the functions F and G, ∀n = 1, . . . , N :

Σ (q+1) n (s (q+1) n ) = 1 M m=1 α (q) mn M m=1 α (q) mn (f m -F(s (q+1) n ))(f m -F(s (q+1) n )) ⊤ , (31)
Γ (q+1) n (s (q+1) n ) = 1 K k=1 β (q) kn K k=1 β (q) kn (g k -G(s (q+1) n ))(g k -G(s (q+1) n )) ⊤ . ( 32
)
For every n = 1, . . . , N , s (q+1) n

is the parameter vector such that:

s (q+1) n = argmax s Q (q) n (s), (33)
where

Q (q) n (s) = - M m=1 α (q) mn ( f m -F(s) 2 Σn(s) + log |Σ n (s)|) - - K k=1 β (q) kn ( g k -G(s) 2 Γn(s) + log |Γ n (s)|). ( 34
)
We stress that the covariances Σ n (s) and Γ n (s) in ( 31) and ( 32) are considered as functions of s ∈ S. Hence, at each iteration of the algorithm, the overall update of the tying parameters can be split into N identical optimization tasks of the form (34). These tasks can be solved in parallel. In general, F and G are non-linear transformations and hence there is no simple closed-form expression for the estimation of the tying parameters.

Generalized EM for Conjugate Mixture Models

The initial parameters selection of the proposed EM algorithm for conjugate mixture models uses the procedure Initialize that is given in Section 6. The maximization step uses two procedures, referred to as Choose and Local Search which are explained in detail in Sections 4 and 5. To determine the number of objects we define the procedure Select that is derived in Section 7. The overall EM procedure is outlined below:

1. Apply procedure Initialize to initialize the parameter vector:

θ (0) = {π (0) 1 , . . . , π (0)
N +1 , λ

1 , . . . , λ

N +1 , s (0)
1 , . . . , s

N , Σ

1 , . . . , Σ

N , Γ

1 , . . . , Γ

N };

2. E step: compute Q(θ, θ (q) ) using equations ( 25) to (28);

3. M step: estimate θ (q+1) using the following sub-steps:

(a) The priors. Compute π

(q+1) 1 , . . . , π (q+1)
N +1 and λ

(q+1) 1 , . . . , λ (q+1)
N +1 using ( 29) and ( 30);

(b) The tying parameters. For each n = 1 . . . N :

• Apply procedure Choose to determine an initial value, denoted by s(0)

n , as proposed in Section 5;

• Apply procedure Local Search to each Q (q) n (s) as defined in (34) start- ing from s(0)
n and set the result to s (q+1) n using the eq. ( 35) specified below;

(c) The covariance matrices. For every n = 1 . . . N , use ( 31) and ( 32 • Initialize: this procedure aims at providing the initial parameter values θ (0) . Its performance has a strong impact on the time required for the algorithm to converge. In Section 6 we propose different initialization strategies based on singlespace cluster detection.

• Select: this procedure applies the BIC-like criterion to determine the number of objects N . In Section 7 propose the consistent criterion for the case of conjugate mixture models.

• Choose: the goal of this procedure is to provide at each M step initial values s(0) 1 , . . . , s(0) N which are likely to be close to the global maxima of the functions 34). The exact form of this procedure is important to ensure the ability of the subsequent Local Search procedure to find these global maxima. We will use results on global search algorithms [START_REF] Zhigljavsky | Stochastic Global Optimization[END_REF] and propose different variants in Section 5.

Q (q) n (s) in (
• Local Search: an important requirement of this procedure is that it finds a local maximum of the

Q (q)
n (s)'s starting from any arbitrary point in S. In this work, we will consider procedures that consist in iterating a local update of the form (ν is the iteration index):

sn (ν+1) = sn (ν) + H (q,ν) n ∇Q (q) n ( sn (ν) ), (35)
with H (q,ν) n being a positive definite matrix that may vary with ν. When the gradient ∇Q (q) n (s) is Lipschitz continuous with some constant L (q) n , an appropriate choice that guarantees the increase of

Q (q) n (s (ν) ) at each iteration ν, is to choose H (q,ν) n such that it verifies H (q,ν) n ≤ 2/L (q)
n . Different choices for H (q,ν) n are possible and they correspond to different optimization methods that belong, in general, to the variable metric class. For example H (q,ν) n = 2 L (q) n I leads to gradient ascent, while taking H (q,ν) n as a scaled inverse of the Hessian matrix would lead to a Newton-Raphson optimization step. Other possibilities include Levenberg-Marquardt and quasi-Newton methods.

Analysis of the Local Search Procedure

Each instance of (34) for n = 1, . . . , N can be solved independently. In this section we focus on providing a set of conditions under which each iteration of our algorithm guarantees that the objective function 34) is increased. We start by rewriting (34) more conveniently in order to perform the optimization with respect to s ∈ S. To simplify the notation, the iteration index q is sometimes omitted. We simply write

Q (q) n (s) in (
Q n (s) for Q (q) n (s). RR n°7117 Let ᾱn = M m=1 α (q) mn and βn = K k=1 β (q)
kn denote the average object weights in each one of the two modalities. We introduce α n = ᾱ-1 n (α

(q) 1n , . . . , α (q)
Mn ) and

β n = β-1 n (β (q) 1n , . . . , β (q)
Kn ) the discrete probability distributions obtained by normalizing the object weights. We denote by F and G the random variables that take their values in the discrete sets {f 1 , . . . , f m , . . . , f M } and {g 1 , . . . , g k , . . . , g K }. It follows that the expressions for the optimal variances (31) and ( 32) as functions of s, can be rewritten as:

Σ (q+1) n (s) =E αn [ F -F(s) F -F(s) ⊤ ], (36)
Γ (q+1) n (s) =E βn [ G -G(s) G -G(s) ⊤ ], (37)
where E αn and E βn denote the expectations with respect to the distributions α n and β n . Using some standard projection formula, it follows that the covariances are:

Σ (q+1) n (s) =V f + v f v ⊤ f , (38)
Γ (q+1) (s) =V g + v g v ⊤ g , (39)
where V f and V g are the covariance matrices of F and G respectively under distributions α n and β n , and v f and v g are vectors defined by:

v f = E αn [F ] -F(s), (40) v g = E βn [G] -G(s). (41)
For convenience we omit the index

n for V f , V g , v f and v g . Let f n = E αn [F ] and ḡn = E βn [G]
. This yields:

f n = ᾱ-1 n M m=1 α (q) mn f m , (42)
ḡn = β-1 n K k=1 β (q) kn g k , (43)
V f = ᾱ-1 n M m=1 α (q) mn f m f ⊤ m -f n f ⊤ n , (44)
V g = β-1 n K k=1 β (q) kn g k g ⊤ k -ḡn ḡ⊤ n . ( 45
)
Next we derive a simplified expression for Q n (s) in ( 34) in order to investigate its properties. Notice that one can write (34) as the sum

Q n (s) = Q n,F (s) + Q n,G (s), with: Q n,F (s) = - M m=1 α (q) mn ( f m -F(s) 2 Σ (q+1) n (s) + log |Σ (q+1) n (s)|), (46)
INRIA and a similar expression for Q n,G (s). Eq. ( 46) can be written:

Q n,F (s) = -ᾱn (E αn [(F -F(s)) ⊤ Σ (q+1) n (s) -1 (F -F(s))] + log |Σ (q+1) n (s)|).
(47) The first term of (47) can be further divided into two terms:

E αn [(F -F(s)) ⊤ Σ (q+1) n (s) -1 (F -F(s))] = =E αn [(F -f n ) ⊤ Σ (q+1) n (s) -1 (F -f n )] + v ⊤ f Σ (q+1) n (s) -1 v f . ( 48
)
The Sherman-Morrison formula applied to (38) leads to

Σ (q+1) n (s) -1 = V -1 f -V -1 f v f v ⊤ f V -1 f /(1 + D n,F (s)), (49)
with:

D n,F (s) = F (s) -f n 2 V f . ( 50
)
It follows that ( 48) can be written as the sum of:

E αn [(F -f n ) ⊤ Σ (q+1) n (s) -1 (F -f n )] = C f - D n,F (s) 1 + D n,F (s) , (51)
and of

v ⊤ f Σ (q+1) n (s) -1 v f = D n,F (s) 1 + D n,F (s) . (52)
Hence the first term of (47), namely ( 48) is equal to C f which is constant with respect to s. Moreover, applying the matrix determinant lemma to the second term of (47) we successively obtain:

log |Σ (q+1) n (s)| = log |V f + v f v ⊤ f | = log |V f | + log(1 + v ⊤ f V -1 f v f ) = = log |V f | + log(1 + D n,F (s)). ( 53
)
It follows that there is only one term depending on s in (47):

Q n,F (s) = -ᾱn (C f + log |V f | + log(1 + D n,F (s))) . (54)
Repeating the same derivation for the second sensorial modality we obtain the following equivalent form of (34):

Q n (s) = -ᾱn log(1 + D n,F (s)) -βn log(1 + D n,G (s)) + C, ( 55
)
where C is some constant not depending on s.

Using this form of Q n (s), we can now investigate the properties of its gradient ∇Q n (s). It appears that under some regularity assumptions on F and G, the gradient ∇Q n (s) is bounded and Lipschitz continuous. The corresponding theorem is formulated and proved. First we establish as a lemma some technical results, required to prove the theorem. In what follows, for any matrix V, the matrix norm used is the operator norm V = sup v =1 Vv . For simplicity, we further omit the index n.

Lemma 1. Let V be a symmetric positive definite matrix. Then the function

ϕ(v) = Vv /(1 + v ⊤ Vv) is bounded by ϕ(v) ≤ C ϕ (V) with C ϕ (V) = V /2 and is Lipschitz continuous: ∀v, ṽ ϕ(v) -ϕ(ṽ) ≤ L ϕ (V) v -ṽ , where L ϕ (V) = V (1 + µ(V)/2) is the Lipschitz constant and µ(V) = V V -1
is the condition number of V.

Proof: We start by introducing w = Vv so that ϕ(v) = φ(w) = w /(1 + w ⊤ V -1 w). As soon as w ⊤ V -1 w ≥ λ min w 2 (where we denoted by λ min the smallest eigenvalue of V -1 , so that in fact λ min = V -1 ), to find the maximum of φ(w) we should maximize the expression t/(1 + λ min t 2 ) for t = w ≥ 0. It is reached at the point t * = λ -1/2 min . Substituting this value into the original expressions gives

ϕ(v) ≤ V /2.
To compute the Lipschitz constant L ϕ we consider the derivative:

∇ φ′ (w) = (1 + w ⊤ V -1 w)w -2 w 2 V -1 w w (1 + w ⊤ V -1 w) 2 ≤ 1 + 2 V -1 w 2 (1 + w ⊤ V -1 w) 2 ,
from where we find that ∇ φ′ (w) ≤ 1+µ(V)/2, and so

L ϕ = V (1+µ(V)/2).
This lemma yields the following main result for the gradient ∇Q: Theorem 1. Assume functions F and G and their derivatives F ′ and G ′ are Lipschitz continuous with constants L F , L G , L ′ F and L ′ G respectively. Then the gradient ∇Q is bounded and Lipschitz continuous with some constant L.

Proof: From (55) the gradient ∇Q can be written as:

∇Q(s) = ∇Q F (s) + ∇Q G (s) = = 2 ᾱF ′ ⊤ (s)V -1 f ( f -F(s)) 1 + D F (s) + 2 βG ′ ⊤ (s)V -1 g (ḡ -G(s)) 1 + D G (s) . (56) It follows from Lemma 1 that ∇Q F (s) ≤ 2L F ᾱC ϕ (V -1 f ) and ∇Q G (s) ≤ 2L G βC ϕ (V -1 g ).
The norm of the gradient is then bounded by:

∇Q(s) ≤ 2L F ᾱC ϕ (V -1 f ) + 2L G βC ϕ (V -1 g ). (57)
Considering the norm ∇Q F (s) -∇Q F (s) , we introduce v 1 = f -F(s) and v 2 = f -F(s). Then we have:

∇Q F (s) -∇Q F (s) ≤ 2 ᾱ (F ′ (s) -F ′ (s)) ⊤ V -1 f v 1 1 + v 1 2 V f + + F ′ ⊤ (s)V -1 f v 2 1 + v 2 2 V f - F ′ ⊤ (s)V -1 f v 1 1 + v 1 2 V f . ( 58
)

INRIA

Using Lemma 1 with V -1 f we have:

∇Q F (s) -∇Q F (s) ≤ 2 ᾱ L ′ F C ϕ (V -1 f ) + L 2 F L ϕ (V -1 f ) s -s .
The same derivations can be performed for ∇Q G (s), so that finally we get:

∇Q G (s) -∇Q G (s) ≤ L s -s , (59)
where the Lipschitz constant is given by:

L = 2 ᾱ L ′ F C ϕ (V -1 f ) + L 2 F L ϕ (V -1 f ) + 2 β L ′ G C ϕ (V -1 g ) + L 2 G L ϕ (V -1 g ) . ( 60
)
To actually construct the non-decreasing sequence in (35), we make use of the following fundamental result on variable metric gradient ascent algorithms.

Theorem 2 [START_REF] Polyak | Introduction to Optimization[END_REF]). Let the function Q : R d → R be differentiable on R d and its gradient ∇Q be Lipschitz continuous with constant L. Let the matrix H be positive definite, such that H ≤ 2 L . Then the sequence Q(s (ν) ), defined by s(ν+1) = s(ν) + H∇Q(s (ν) ) is non-decreasing.

This result shows that for any functions F and G that verify the conditions of Theorem 1, using ( 35) with H = 2 L I, we are able to construct a non-decreasing sequence and an appropriate Local Search procedure. Notice however, that its guaranteed theoretical convergence speed is linear. It can be improved in several ways.

First, the optimization direction can be adjusted. For certain problems, the matrix H can be chosen as in variable metric algorithms, such as Newton-Raphson method, quasi-Newton methods or Levenberg-Marquardt method, provided that it satisfies the conditions of Theorem 2. Second, the optimization step size can be increased based on local properties of the target function. For example, at iteration ν, if when considering the functions F and G on some restricted domain S (ν) there exist smaller local Lipschitz constants L (ν)

F , L (ν) G , L ′(ν) F and L ′(ν) G , H can be set to H = 2 L (ν) I with L (ν) smaller than L. It follows that s(ν+1) -s(ν) ≤ 2 L (ν) ∇Q(s (ν) )
, which means that one can take the local constants, L (ν)

F , L (ν) G , L ′(ν) F and L ′(ν) G if they are valid in the ball B ρ (ν) (s (ν) ) with ρ (ν) = 2 L (ν) 2L (ν) F ᾱC ϕ (V -1 f ) + 2L (ν) G βC ϕ (V -1 g ) . ( 61
)
5 Global Search and the Choose Procedure

Theorem 1 allows us to use the improved global random search techniques for Lipschitz continuous functions [START_REF] Zhigljavsky | Theory of Global Random Search[END_REF]. These algorithms are known to converge, in the sense that generated point sequences fall infinitely often into an arbitrarily small neighbourhood of the optimal points set. For more details and convergence conditions see Theorem 3.2.1 and the discussion that follows in [START_REF] Zhigljavsky | Theory of Global Random Search[END_REF]. A proper choice of the initial value s(0) not only guarantees to find the global maximum, but can also be used to increase the convergence speed. A basic strategy is to draw samples in S, according to some sequence of distributions over S, that verifies the convergence conditions of global random search methods. However, the speed of convergence of such an algorithm is quite low.

Global random search methods can also be significantly improved by taking into account some specificities of the target function. Indeed, in our case, function ( 55) is made of two parts for which the optimal points are known and are respectively f and ḡ. If there exists s(0) such that s(0) ∈ F -1 ( f ) ∩ G -1 (ḡ), then it is the global maximum and the M step solution is found. Otherwise, one can sample S in the vicinity of the set F -1 ( f ) ∪ G -1 (ḡ) to focus on a subspace that is likely to contain the global maximum. This set is, generally speaking, a union of two manifolds. For sampling methods on manifolds we refer to [START_REF] Zhigljavsky | Theory of Global Random Search[END_REF]. An illustration of this technique is given in Section 8.

Another possibility is to use a heuristic that function (55) does not change much after one iteration of the EM algorithm. Then, the initial point s(0) for the current iteration can be set to the optimal value computed at the previous iteration. However, in general, this simple strategy does not yield the global maximum, as can be seen from the results in Section 9.1.

Algorithm Initialization and the Initialize Procedure

In this section we focus on the problem of selecting the initial values θ (0) for the model parameters. As it is often the case with iterative optimization algorithms, the closer θ (0) is to the optimal parameter values, the less time the algorithm would require to converge. Within the framework of conjugate mixture models we formulate two main strategies, namely Observation Space Candidates and Parameter Space Candidates that attempt to find a good initialization.

The Observation Space Candidates strategy consists in searching for cluster centers in single modality spaces F and G to further map them into the parameter space S, and select the best candidates. More specifically, we randomly select an observation f m (or g k ) and run the mean shift algorithm [START_REF] Comaniciu | Mean shift: a robust approach towards feature space analysis[END_REF] in the corresponding space to find local modes of the distribution, which are called candidates. The sets of candidate points { fi } i∈I and {ĝ j } j∈J are further rarefied, that is if f i1 -f i2 ≤ ε f for some i 1 = i 2 and for some threshold ε > 0, we eliminate one of these points. These rarefied sets are then mapped to S. If one of the observation space mappings, for example F, is non-injective, for each f i we need to select a point s i ∈ F -1 ( f i ) that is the best in some sense. We consider observations density in the other observation spaces around an image of s i as the optimality measure of s i . This can be estimated through calculation of the k-th nearest neighbour distance (k-NN) in INRIA the corresponding observation space. The final step is to choose N points out of these candidates to initialize the cluster centers {s 1 , . . . , s N }, so that the inter-cluster distances are maximized. This can be done using, for example, hierarchical clustering. The variances Σ 1 , . . . , Σ N and Γ 1 , . . . , Γ N are then calculated by standard empirical variance formulas based on observations, that are closest to the corresponding class center. The priors π 1 , . . . , π N +1 and λ 1 , . . . , λ N +1 are set to be equal.

The Parameter Space Candidates strategy consists in mapping all the observations to the parameter space S, and performing subsequent clustering in that space. More specifically, for every observation f m and g k we find an optimal point from the corresponding preimage F -1 (f m ) and G -1 (g k ). The optimality condition is the same as in the previous strategy, that is we compare the local observation densities using k-NN distances. Then one proceeds with selecting local modes in space S using the meanshift algorithm, and initializing N cluster centers {s 1 , . . . , s N } from all the candidates thus calculated. The estimation of variances and priors is exactly the same as in the previous strategy.

The second strategy proved to be better when performing simulations (see Section 9). This can be explained by possible errors in finding the preimage of an observation space point in the parameter space. Thus mapping a rarefied set of candidates to the parameter space is less likely to make a good guess in that space than mapping all the observations and finding the candidates directly in the parameter space.

Estimating the Number of Components and the Select Procedure

To choose the N that best corresponds to the data, we perform model selection based on a criterion that resembles the BIC criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF]. We consider the score function of the form

BIC N = -2L(f, g, θN ) + D N log(M + K), ( 62
)
where θN is the ML estimate obtained by the proposed EM algorithm, L(f, g, θ) is given by ( 15) and

D N = N d + 2 + 1 2 (r 2 + p 2 + r + p) is the dimensionality of the model.
As in the case of (non-conjugate) Gaussian mixture models, we cannot derive the criterion from the Laplace approximation of the probability P (f, g|N = N 0 ) because of the Hessian matrix of L(f, g, θ) that is not necessarily positive definite [START_REF] Aitkin | Estimation and hypothesis testing in finite mixture models[END_REF][START_REF] Quinn | A note on the Aitkin-Rubin approach to hypothesis testing in mixture models[END_REF]. Nevertheless, we can use the same arguments as those used in [START_REF] Keribin | Consistent estimation of the order of mixture models[END_REF] for Gaussian mixture models to show that the criterion is consistent, i.e. if N * is the number of components in the real model that generated f and g, then

N BIC → N * a.s., when M, K → ∞, (63)
provided variances Σ 1 , . . . , Σ N , Γ 1 , . . . , Γ N are non-degenerate and the sequence M M+K has only one accumulation point (i.e. has a limit).

The BIC-like criterion (62) shows good performance on both simulated and real data (see Section 9), choosing correctly the number of objects in all the cases.

Clustering Using Auditory and Visual Data

We illustrate the method in the case of audiovisual (AV) objects. Objects could be characterized both by their locations in space and by their auditory status, i.e., whether they are emitting sounds or not. These object characteristics are not directly observable and hence they need to be inferred from sensor data, e.g., cameras and microphones. These sensors are based on different physical principles, they operate with different bandwidths and sampling rates, and they provide different types of information. On one side, light waves convey useful visual information only indirectly, on the premise that they reflect onto the objects' surfaces. A natural scene is composed of many objects/surfaces and hence the task of associating visual data with objects is a difficult one. On the other side, acoustic waves convey auditory information directly from the emitter to the receiver but the observed data is perturbed by the presence of reverberations, of other sound sources, and of background noise. Moreover, very different methods are used to extract information from these two sensor types. A wide variety of computer vision principles exist for extracting 3D points from a single image or from a pair of stereoscopic cameras [START_REF] Forsyth | Computer Vision -A Modern Approach[END_REF] but practical methods are strongly dependent on the ligthing conditions and on the properties of the objects' surfaces (presence or absence of texture, color, shape, reflectance, etc.). Similarly, various algorithms were developed to locate sound sources using a microphone pair based on interaural time differences (ITD) and on interaural level differences (ILD) (Wang and Brown, 2006;[START_REF] Christensen | Integrating pitch and localisation cues at a speech fragment level[END_REF], but these cues are difficult to interpret in natural settings due to the presence of background noise and of other reverberant objects. A notable improvement consists in the use a larger number of microphones [START_REF] Dibiase | Robust localization in reverberant rooms[END_REF]. Nevertheless, the extraction of 3D sound source positions from several microphone observations results in inaccurate estimates. We show below that our method can be used to combine visual and auditory observations to detect and localize objects. A typical example where the conjugate mixture models framework may help is the task of locating several speaking persons.

Using the same notations as above, we consider two sensor spaces. The multimodal data consists of M visual observations f and of K auditory observations g. We consider data that are recorded over a short time interval [t 1 , t 2 ], such that one can reasonably assume that the AV objects have a stationary spatial location. Nevertheless, it is not assumed here that the AV objects, e.g., speakers, are static: lip movements, head and hand gestures are tolerated. We address the problem of estimating the spatial locations of all the objects that are both seen and heard. Let N be the number of objects and in this case each object is described by a three dimensional parameter vector

s n = (x n , y n , z n ) ⊤ .
The AV data are gathered using a pair of stereoscopic cameras and a pair of omnidirectional microphones, i.e., binocular vision and binaural hearing. A visual observation INRIA vector f m = (u m , v m , d m ) ⊤ corresponds to a 2D image location (u m , v m ) and to an associated binocular disparity d m . Considering a projective camera model [START_REF] Faugeras | Three Dimensional Computer Vision: A Geometric Viewpoint[END_REF] it is straightforward to define an invertible function F : R 3 → R 3 that maps s = (x, y, z) ⊤ onto f = (u, v, d) ⊤ :

F(s) = x z , y z , 1 z ⊤ and F -1 (f ) = u d , v d , 1 d ⊤ . ( 64
)
This model corresponds to a rectified camera pair [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] and it can be easily generalized to more complex binocular geometries (Hansard andHoraud, 2008, 2007). Without loss of generality one can use a sensor-centered coordinate system to represent the object locations.

Similarly one can use the auditory equivalent of disparity, namely the interaural time difference (ITD) widely used by auditory scene analysis methods (Wang and Brown, 2006). The function G : R 3 → R maps s = (x, y, z) ⊤ onto a 1D audio observation:

g = G(s) = 1 c s -s M1 -s -s M2 . ( 65
)
Here c is the sound speed and s M1 and s M2 are the 3D locations of the two microphones in the sensor-centered coordinate system. Each isosurface defined by ( 65) is represented by one sheet of a two-sheet hyperboloid in 3D. Hence, each audio observation g constrains the location of the auditory source to lie onto a 2D manifold.

In order to perform audiovisual clustering based on the conjugate EM algorithm, Theorem 1 (Section 4) must hold for both ( 64) and ( 65), namely the functions F and G and their derivatives are Lipschitz continuous. We prove the following theorem: Theorem 3. The functions F, F ′ , G and G ′ are Lipschitz continuous with constants

L F = z -1 min √ 3, L ′ F = z -2 min , L G = s M1 -s M2 (cR) -1 and L ′ G = 3(cR) -1 in the domain S = {|z| > z min > 1} ∩ min{ s -s M1 , s -s M2 } > R > 1 .
Proof: The derivatives of F and G are given by:

F ′ (s) = 1 z   1 0 -x/z 0 1 -y/z 0 0 -1/z   (66) G ′ (s) = 1 c s -s M1 s -s M1 - s -s M2 s -s M2 . ( 67
)
The eigenvalues of F ′ (s) are 1/z and -1/z 2 , so

F ′ (s) ≤ max{z -1 , z -2 } ≤ z -1
min , from which it follows that L F can be taken as

L F = z -1 min √ 3. Also F ′ (s) - F ′ (s) ≤ max{|z -1 -z-1 |, |z -2 -z-2 |} ≤ z -2 min s -s , so that L ′ F can be set to L ′ F = z -2 min . Introducing e 1 = s-sM1 s-sM 1 and e 2 = s-sM2 s-sM 2 , it comes e 1 = e 2 = 1
and G ′ (s) = 1 c (e 1e 2 ). Provided that ss M1 and ss M2 are both greater than R, it follows G ′ (s) = 1 c e 1e 2 ≤ s M1s M2 (cR) -1 and so L G = s M1s M2 (cR) -1 . Then, the second derivative of G is given by

G ′′ (s) = 1 c s -s M1 (I -e 1 e ⊤ 1 ) - 1 c s -s M2 (I -e 2 e ⊤ 2 )
.

so that G ′′ (s) ≤ 1 c s-sM1 - 1 c s-sM2 + sup v =1 2e1e ⊤ 1 v c min{ s-sM1 , s-sM2 } ≤ 3(cR) -1 ,
and L ′ G can be set to L ′ G = 3(cR) -1 . This result shows that under some natural conditions (The AV objects should not be too close to the sensors) the conjugate EM algorithm described in Section 3.3 can be applied. The constant L given by Lemma 1 guarantees a certain (worst-case) convergence speed. In practice, we can use the techniques mentioned in Sections 4 and 5 to accelerate the algorithm. First, to speed up the local optimization step, local Lipschitz constants can be computed based on the current value of parameter s(ν) . Equation ( 61) gives the largest possible step size ρ (ν) , so setting z ν) , provides local Lipschitz constants that insure the update not to quit ν) . Second, we propose four possibilities to set the initial object parameter values s(0) n : (i) it can be taken to be the previously estimated object position s (q-1) n , (ii) it can be set to F -1 ( f ) (as soon as F is injective in S), (iii) it can be found through sampling of the manifold G -1 (ḡ) by selecting the sampled value which gives the largest Q value, or (iv) similarly through sampling directly in S. Comparisons are reported in the following sections.

(ν) min = z (ν) -ρ (ν) and R (ν) = min{ s(ν) -s M2 , s(ν) -s M1 }-ρ (
S (ν) = {|z| > z (ν) min } ∩ min{ s -s M1 , s -s M2 } > R (

Experimental Validation

Experiments with Simulated Data

Our algorithm is first illustrated on simulated data. For simplicity we consider (u, d) and (x, z) coordinates so that F ⊆ R 2 and S ⊆ R 2 . Notice however that this preserves the projective nature of the mapping F, it does not qualitatively affect the results and allows to better understand the algorithm performance. We consider three objects defined in S by s n , n = 1, 2, 3. We simulated three cases: well-separated objects (GoodSep), partially occluded objects (PoorSep) and poor precision in visual observations for well-separated objects (PoorPrec). The ground-truth object locations (x, z) for the GoodSep and PoorPrec cases are the same, namely s 1 = (-300, 1000), s 2 = (10, 800) and s 3 = (500, 1500). In the PoorSep case, the coordinates are respectively s 1 = (-300, 1000), s 2 = (10, 800) and s 3 = (100, 1500). The data in both observation spaces F and G was simulated from a mixture model with three Gaussian components and a uniform component that models the outliers. The means of the Gaussian components are computed using F(s n ) and G(s n ), n = 1, 2, 3. An example of simulated data for the three mentioned configurations is shown in Figure 2 Initialization. We compared two strategies, Observation Space Candidates (OSC) and Parameter Space Candidates (PSC) that are proposed in Section 6. Their performance is summarized in Figure 3. It shows the mean and variance of the likelihood value L(f, g, θ) for initial parameters θ

OSC and θ

PSC chosen by OSC and PSC strategies respectively. For the total number of clusters N = 1, . . . , 5 and different object configurations, we calculate the statistics based on 10 initializations. The analysis shows that the PSC strategy performs at least as well as the OSC strategy, or even better in some cases. Our explanation is that mappings from observation spaces to parameter space are subject to absolute (and in our case bounded) noise. Mapping all the observations and calculating a candidate point in the parameter space has an averaging effect and reduces the absolute error, compared to the strategy with candidate calculation being performed in an observation space with subsequent mapping to the parameter space. Therefore in what follows, all the results are obtained based on the PSC initialization strategy.

Optimization. We compared several versions of the algorithm based on various Choose and Local Search strategies. For the initial values sn (0) , we considered the following possibilities: the optimal value computed at a previous run of the algorithm (IP), the value predicted from visual data (IV), the value predicted from audio data (IA) and the value obtained by global random search (IG). More specifically: • When initializing from visual data (IV), the average value f n , calculated in the current E-step of the algorithm for every n, was mapped to the parameter space and sn (0) set to sn (0) = F -1 ( f n ) using the injectivity of F.

• When initializing from audio data (IA), G -1 (ḡ n ) defines a manifold. The general strategy here would be to find the optimal point that lies on this surface. We achieved this through random search based on a uniform sampling on the corresponding part of the hyperboloid (see [START_REF] Zhigljavsky | Theory of Global Random Search[END_REF] for details on sampling from an arbitrary distribution on a manifold); in our experiments we used 50 samples to select the one providing the largest Q (likelihood) value.

• The most general initialization scheme (IG) was implemented using global random search in the whole parameter space S; 200 samples were used in this case.

Local optimization was performed either using basic gradient ascent (BA) or the locally accelerated gradient ascent (AA). The latter used the local Lipschitz constants to augment the step size, as described in Section 4.

Each algorithm run consisted of 70 iterations of the EM algorithm with 10 nondecreasing iterations during the M step.

To check the convergence speed of different versions of the algorithm for the three object configurations we compared the likelihood evolution graphs that are presented in Figure 4. Each graph contains several curves that correspond to five different versions of the algorithm. The acronyms we use to refer to the different versions (for example, IPAA) consist of two parts encoding the initialization (IP) and the local optimization (AA) types. The black dashed line on each graph shows the 'ground truth' likelihood level, that is the likelihood value for the parameters used to generate the data. The meaning of the acronyms is recalled in Table 1.

As expected, the simplest version IPBA that uses none of the proposed acceleration techniques appears to be the slowest. The other variants using basic gradient ascent are then not reported. Predicting a single object parameter value from visual observations (IVAA) does not give any improvement over IPAA, where s(0) is taken from the previous EM iteration. When s(0) is obtained by sampling the hyperboloid predicted from audio observations (IAAA), a significant impact on the convergence speed is observed, especially on early stages of the algorithm, where the predicted value can be quite far from the optimal one. However, 'blind' sampling of the whole parameter space does not bring any advantage: it is much less efficient regarding the number of samples required for the same precision. This suggests that in the general case, the best strategy would be to sample the manifolds F -1 ( f n ) and G -1 (ḡ n ) with possible small perturbations to find the best s(0) estimate and to perform an accelerated gradient ascent afterwards (IAAA). We note that IAAA succeeds in all the cases to find parameter values that are well-fitted to the model in terms of likelihood function (likelihood is greater or equal than that of real parameter values).

Parameter evolution trajectories for the IAAA version of the algorithm in the Good-Sep case are shown in Figures 56. The estimate changes are reflected by the node sizes (from smaller to bigger) and colours (from darker to lighter). The final values are very close to the real cluster centers in all three audio, visual and object spaces. The convergence speed is quite dependent on the initialization. In the provided example the algorithm spent almost a half of useful iterations to disentangle the estimates trying to decide which one corresponds to which class. Another possibility here would be to predict the initial values through sampling in the audio domain. We demonstrate this strategy further when working with real data.

We compared the performance of our algorithm for the three object configurations. For each of them, we computed absolute and relative errors for the object parameter estimations in the different coordinate systems (object, audio and visual spaces). The averages were taken over 10 runs of the algorithm for different PSC initializations, as described above. The results are reported in Table 2. We give object location estimates ŝ = (x, ẑ), f = (û, d) and ĝ in parameter, visual and audio spaces respectively. It appears that the localization precision is quite high. In a realistic setting such as that of Section 9.2, the measurement unit can be set to a millimeter. In that case, the observed precision, in a well-separated objects configuration, it is at worse about 6cm. However, precision in the z coordinate is quite sensible to the variance of the visual data and the object configuration. To get a better idea of the relationship between the variance in object space and the variance in visual space, F -1 can be replaced by its linear approximation given by a first order Taylor expansion. Assuming then that visual data are distributed according to some probability distribution with mean µ F and variance Σ F , it follows that through the linear approximation of F -1 , the variance in object

space is ∂F -1 (µF ) ∂f Σ F ∂F -1 (µF ) ∂f ⊤
. Then, the z coordinate covariance for an object n is approximately proportional to the d covariance for the object multiplied by z 4 n . For distant objects, a very high precision in d is needed to get a satisfactory precision in z. At the same time we observe that the likelihood of the estimate configuration often exceeds the likelihood for real parameter values. This suggests that the model performs well for the given data, but cannot get better precision than that imposed by the data. Selection. To select the optimal number of clusters N we applied the BIC criterion ( 62

f = (u, d) f = (û, d) e a = f -f e r = f -f / f GoodSep Object 1 (-0.

Experiments with Real Data

We evaluated the ability of our algorithms to estimate the 3D locations of AV objects in a person localization task. We considered a typical 'meeting' situation. The data was taken from a database of realistic AV scenarios described in detail in [START_REF] Arnaud | The CAVA corpus: Synchronized stereoscopic and binaural datasets with head movements[END_REF]. A mannequin, with a pair of microphones fixed into its ears and a pair of stereoscopic cameras mounted onto its forehead, served as the acquisition device (the setup was developed and constructed within the POP1 project). The reason for choosing this configuration was to record data from the perspective of a person, i.e. to try to capture what a person would both hear and see while being in a natural AV environment. Each recorded scenario comprised two audio tracks and two image sequences, together with the sensor calibration information. The data we use here is a meeting scenario2 , shown in Figure 8. There are five persons sitting around a table, but only three persons are visible. The recording lasts 25 seconds and contains a total of about 8000 visual and 600 audio observations.

Audio and visual observations were collected over the recording using the following techniques. A standard procedure was used to identify 'interest points' in the left and right images [START_REF] Harris | A combined corner and edge detector[END_REF]. These features were put into binocular correspondence by comparing the local image-structure at each of the candidate points, as described in [START_REF] Hansard | Patterns of binocular disparity for a fixating observer[END_REF]. The cameras were calibrated using Intel's OpenCV3 in order to define the (u, v, d) ⊤ to (x, y, z) ⊤ mapping (64). Auditory disparities were obtained through the analysis of cross-correlogram of the filtered left and right microphone signals for every frequency band [START_REF] Christensen | Integrating pitch and localisation cues at a speech fragment level[END_REF]. To be able to introduce the common parameter space S, we performed audio-visual calibration which consisted in finding microphone positions in camera-related 3D frame.

The meeting situation corresponds to a well-separated targets case (GoodSep simulation). The likelihood evolution reported in Figure 9 shows that convergence is reached in about 20 iterations which is three times faster than in the simulated GoodSep case of Figure 4. This suggests that the proposed initialization strategy is appropriate. Also the 3D position estimates are quite accurate, in particular the natural alignment of the speakers along the table is clearly seen in the XZ domain of the object space. Even though the outlier class is not uniform and the clusters themselves are not Gaussian, our model performed quite well, which illustrates robustness of the proposed model to actual observation-space data distributions.

Conclusions

We proposed a novel framework, conjugate mixture models to cluster heterogeneous data gathered with physically different sensors. Our approach differs from other existing approaches in that it combines in a single statistical model a number of clustering tasks while ensuring the consistency of their results. In addition, the fact that the clustering is performed in observation spaces allows one to get useful statistics on the data, which is an advantage of our approach over particle filtering models. The task of simultaneous clustering in spaces of different nature, related through known functional dependencies to a common parameter space, was formulated as a likelihood maximization problem. Using the ideas underlying the classical EM algorithm we built the conjugate EM algorithm to perform the multimodal clustering task, while keeping attractive convergence properties. The analysis of the conjugate EM algorithm and, more specifically, of the optimization task arising in the M-step, revealed several possibilities to increase the convergence speed. We proposed to decompose the M-step into two procedures, namely the Local Search and Choose procedures, which allowed us to derive a number of acceleration strategies. We exhibited appealing properties of the target function which induced several implementations of these procedures resulting in a significantly improved convergence speed. We introduced the Initialize and Select procedures to efficiently choose initial parameter values and determine the number of clusters in a consistent manner respectively. A non trivial audio-visual localization task was considered to illustrate the conjugate EM performance on both simulated and real data. Simulated data experiments allowed us to assess the average method behaviour in various configurations. They showed that the obtained clustering results were precise as regards the observation spaces under consideration. They also illustrated the theoretical dependency between the precisions in observation and parameter spaces. Real data experiments then showed that the observed data precision was high enough to guarantee high precision in the parameter space.

One of the strong points of the formulated model is that it is open to different useful extensions. It can be easily extended to an arbitrary number J of observation spaces F 1 , . . . , F J . The main results, including Local Search and Choose acceleration strategies stay valid with minor changes. The sum of two terms, related to spaces F and G, would have to be replaced by a sum of J terms corresponding to F 1 , . . . , F J in the formulas of Section 3. Also, the assumption that assignment variables a and b INRIA are independent could be relaxed. An appropriate approach to perform inference in a non independent case would be to consider variational approximations [START_REF] Jordan | An introduction to variational methods for graphical models[END_REF] and in particular a variational EM (VEM) framework. The general idea would be to approximate the joint distribution P (a) by a distribution from a restricted class of probability distributions that factorize as P (a) = M m=1 P (a m ). For any such distribution, our model would be applicable without any changes so that for a variational version of the conjugate EM algorithm, all the results from Section 3 would hold.

It appears that as a generalization of Gaussian mixture models, our model has larger modeling capabilities. It is entirely based on a mathematical framework in which each step is theoretically well-founded. Its ability to provide good results in a non trivial multimodal clustering task is particularly promising for applications requiring the integration of several heterogenious information sources. Therefore, it has advantages over other methods that include ad-hoc processing while being open to incorporation of more task dependent information.

Figure 1 :

 1 Figure 1: Graphical representation of the conjugate mixture model. Circles denote random variables, plates (rectangles) around them represent multiple similar nodes, their number being given in the plates.

  Check for convergence: Terminate, otherwise go to Step 2; 5. Apply procedure Select, use (62) specified below to determine the best N ; INRIA This algorithm uses the following procedures:

Figure 2 :

 2 Figure 2: Simulated data in visual (top) and audio (bottom) observation spaces for three cases: (a) well-separated objects, (b) partially occluded objects, and (c) poor precision of visual observations. The small squares correspond to the ground-truth parameter values. Each one of the two mixtures models (associated with each sensorial modality) contains four components: three objects and one outlier class.

Figure 3 :

 3 Figure 3: Means and variances of log-likelihood values L(f, g, θ) for initial parameters θ (0) OSC and θ (0) PSC chosen by Observation Space Candidates (OSC, red) and Parameter Space Candidates (PSC, blue) strategies respectively, for different numbers of clusters N and different data configurations.

Figure 4 :

 4 Figure 4: Likelihood function evolution for five variants of the algorithm in three cases. Top-left: well-separated objects; top-right: poorly separated objects; bottom: wellseparated object but poor observation precision.

Figure 5 :Figure 6 :

 56 Figure5: IAAA algorithm: parameter evolution and assignment results for the Good-Sep case in audio and visual spaces (note the scale change which corresponds to a zoom on the cluster centers). The initialization (white stars) is based on the PSC strategy. Ground truth means are marked with squares. The evolution is shown by circles from smaller to bigger, from darker to brighter. Observatio ns assignments are depicted by different markers (•, * and × for the three object classes) in visual space and are colour-coded in audio space. Due to the zoom, outliers are not visible on these figures.

Figure 7 :

 7 Figure 7: BIC score graphs for the three object configurations, evaluated for models trained for different total number of clusters N .

Figure 8 :

 8 Figure 8: Visual observations on the left and right camera images. Circles depict the 'interest points', squares show those of them that are matched to some point from the other image. The epipolar lines correspond to a point marked by a star in the opposite image.

Table 1 :

 1 Acronyms used for five variants of the conjugate EM algorithm. Variants correspond to different choices for the Choose and Local search procedures.

	Acronym	s(0) initialization (Choose)	Local optimization (Search)
	IPBA	previous iteration value	basic gradient ascent
	IGAA	global random search	accelerated gradient ascent
	IVAA	predicted value from visual data	accelerated gradient ascent
	IPAA	previous iteration value	accelerated gradient ascent
	IAAA	audio predicted manifold sampling accelerated gradient ascent

Table 2 :

 2 IAAA algorithm: object location estimates in parameter, visual and audio spaces for GoodSep, PoorSep and PoorPrec object configurations. The estimates are calculated based on ten runs of the algorithm with PSC initializations.

		Ground Truth Estimates Mean Absolute Error	Relative Error
	Parameter Space GoodSep Object 1 Object 2 Object 3 PoorSep Object 1 Object 2 Object 3 PoorPrec Object 1 Object 2 Object 3 Visual Space	s = (x, z) (-300, 1000) (10, 800) (500, 1500) (-300, 1000) (10, 800) (100, 1500) (-300, 1000) (10, 800) (500, 1500)	ŝ = (x, ẑ) (-300.13, 997.81) (9.28, 804.46) (513.56, 1555.23) (-307.47, 1028.38) (14.19, 895.69) (105.02, 1447.49) (-208.86, 698.51) (8.44, 703.97) (507.65, 1533.8)	e a = ŝ -s 2.2 4.52 56.86 29.35 95.79 52.75 314.97 96.04 34.66	e r = ŝ -s / s 2.1 • 10 -3 5.7 • 10 -3 3.5 • 10 -2 2.8 • 10 -2 1.2 • 10 -1 3.5 • 10 -2 0.3 1.2 • 10 -1 2.2 • 10 -2

RR n°7117

http://perception.inrialpes.fr/POP/

http://perception.inrialpes.fr/CAVA Dataset/Site/data.html#M1

http://www.intel.com/technology/computing/opencv INRIA

Acknowledgements

The authors would like to thank Miles Hansard (INRIA Grenoble Rhône-Alpes) and Heidi Christensen (Department of Computer Science, University of Sheffield) for providing their software for visual and auditory feature detection.

This work was supported by the Perception-on-Purpose project, under EU grant FP6-IST-2004-027268. * MISTIS team † PERCEPTION team Conjugués pour le Clustering de Données Multimodales

To initialize the parameter values we used the Parameter Space Candidates (PSC) initialization strategy, as described in Section 6. As in the case of simulated data, we used the BIC criterion to select the optimal number of clusters N = 3 for the data. The results are shown in Figure 9.