1

Conjugate Mixture Models for Clustering Multimodal Data111This work was
supported by the Perception-on-Purpose project, under EU grant
FP6-IST-2004-027268.
Neural Computation, Volume 23, Issue 2, February 2011, 517-557

Vasil Khalidov†, Florence Forbes† and Radu Horaud†
†INRIA Grenoble Rhône-Alpes, 655, avenue de l’Europe
38330 Montbonnot Saint-Martin, FRANCE

Keywords: mixture models, EM algorithm, multisensory fusion, audiovisual
integration, Lipschitz continuity, global optimization.

###### Abstract

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 there 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 thoroughly 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.

##  1 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 modelling of
biological multisensory perception. This includes the issues of how a human
detects objects that are both seen and touched (Pouget et al., 2002; Ernst and
Banks, 2002), seen and heard (Anastasio et al., 2000; King, 2004, 2005) 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 (Haykin and
Chen, 2005). More generally, multisensory fusion (Hall and McMullen, 2004;
Mitchell, 2007) is highly relevant in various other research domains, such as
target tracking (Smith and Singh, 2006) based on radar and sonar data (Naus
and van Wijk, 2004; Coiras et al., 2007), mobile robot localization with laser
rangefinders and cameras (Castellanos and Tardos, 1999), robot manipulation
and object recognition using both tactile and visual data (Allen, 1995; Joshi
and Sanderson, 1999), underwater navigation based on active sonar and
underwater cameras (Majumder et al., 2001), audio-visual speaker detection
(Beal et al., 2003; Perez et al., 2004; Fisher III and Darrell, 2004), speech
recognition (Heckmann et al., 2002; Nefian et al., 2002; Shao and Barker,
2008), 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 (Beal et
al., 2003; Kushal et al., 2006; Smith and Singh, 2006). 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
(Castellanos and Tardos, 1999; Majumder et al., 2001). 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 (Fisher III et al., 2001) 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 (Checka et al., 2004; Chen and
Rui, 2004; Gatica-Perez et al., 2007). 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 (Chen and Rui, 2004; Gatica-Perez et al., 2007) 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 in the parameter space, but no statistics are
gathered for the observation spaces. Recently (Hospedales and Vijayakumar,
2008) extended the single-object model of (Beal et al., 2003) 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 (McLachlan and Peel, 2000). We propose a variant of the EM algorithm
(Dempster et al., 1977; McLachlan and Krishnan, 1996) 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 M-step
of the proposed algorithm (Zhigljavsky and Žilinskas, 2008). 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 10 analyses in detail the performances of the
proposed model under various practical conditions with both simulated and real
data. Finally, Section 11 the paper and provides directions for future work.

##  2 Mixture Models for Multimodal Data

We consider $N$ objects $n=1\ldots N$. Each object $n$ is characterized by a
parameter vector of dimension $d$, denoted by
$\mbox{\boldmath$s$}_{n}\in\mathbb{S}\subseteq\mathbb{R}^{d}$. The set
$\mbox{\bf
s}=\\{\mbox{\boldmath$s$}_{1},\ldots,\mbox{\boldmath$s$}_{n},\ldots,\mbox{\boldmath$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
$\mbox{\bf
f}=\\{\mbox{\boldmath$f$}_{1},\ldots,\mbox{\boldmath$f$}_{m},\ldots,\mbox{\boldmath$f$}_{M}\\}$
and $\mbox{\bf
g}=\\{\mbox{\boldmath$g$}_{1},\ldots,\mbox{\boldmath$g$}_{k},\ldots,\mbox{\boldmath$g$}_{K}\\}$
lying in two different observation spaces of dimensions $r$ and $p$,
$\mbox{\boldmath$f$}_{m}\in\mathbb{F}\subseteq\mathbb{R}^{\,r}$ and
$\mbox{\boldmath$g$}_{k}\in\mathbb{G}\subseteq\mathbb{R}^{p}$.

One key ingredient of our approach is that we consider the transformations:

$\left\\{\begin{array}[]{c}\mathcal{F}:\mathbb{S}\to\mathbb{F}\\\
\mathcal{G}:\mathbb{S}\to\mathbb{G}\end{array}\right.$ (1)

that map $\mathbb{S}$ respectively into the observation spaces $\mathbb{F}$
and $\mathbb{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 $\mathcal{F}$ and $\mathcal{G}$ are non-
linear.

An assignment variable is associated with each observation, thus indicating
the object that generated the observation:
$\mbox{\boldmath$A$}=\\{A_{1},\ldots,A_{m},\ldots,A_{M}\\}$ and
$\mbox{\boldmath$B$}=\\{B_{1},\ldots,B_{k},$ $\ldots,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
$\mbox{\boldmath$f$}_{m}$ (resp. $\mbox{\boldmath$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
$\mbox{\boldmath$f$}_{m}$ (resp. $\mbox{\boldmath$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(\mbox{\bf f},\mbox{\bf
g}|\mbox{\boldmath$a$},\mbox{\boldmath$b$})=\prod\limits_{m=1}^{M}P(\mbox{\boldmath$f$}_{m}|a_{m})\prod\limits_{k=1}^{K}P(\mbox{\boldmath$g$}_{k}|b_{k}).$
(2)

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

$P(\mbox{\boldmath$a$},\mbox{\boldmath$b$})=\prod\limits_{m=1}^{M}P(a_{m})\prod\limits_{k=1}^{K}P(b_{k}).$
(3)

As discussed in Section 11, 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 (Jordan et al., 1998; Celeux et al.,
2003).

Next we define the following probability density functions, for all $n=1\ldots
N,N+1$, for all $\mbox{\boldmath$f$}_{m}\in\mathbb{F}$ and for all
$\mbox{\boldmath$g$}_{k}\in\mathbb{G}$:

$\displaystyle P_{n}^{\mathbb{F}}(\mbox{\boldmath$f$}_{m})$ $\displaystyle=$
$\displaystyle P(\mbox{\boldmath$f$}_{m}|A_{m}=n),$ (4) $\displaystyle
P_{n}^{\mathbb{G}}(\mbox{\boldmath$g$}_{k})$ $\displaystyle=$ $\displaystyle
P(\mbox{\boldmath$g$}_{k}|B_{k}=n).$ (5)

More specifically, the likelihoods for an observation to belong to an object
$n$ are Gaussian distributions whose means
$\mathcal{F}(\mbox{\boldmath$s$}_{n})$ and
$\mathcal{G}(\mbox{\boldmath$s$}_{n})$ correspond to the object’s parameter
vector $\mbox{\boldmath$s$}_{n}$ mapped to the observations spaces by the
transformations $\mathcal{F}$ and $\mathcal{G}$:

$\displaystyle P_{n}^{\mathbb{F}}(\mbox{\boldmath$f$}_{m})$ $\displaystyle=$
$\displaystyle{\mathcal{N}}(\mbox{\boldmath$f$}_{m};\;\mathcal{F}(\mbox{\boldmath$s$}_{n}),{\bf\Sigma}_{n}),$
(6) $\displaystyle P_{n}^{\mathbb{G}}(\mbox{\boldmath$g$}_{k})$
$\displaystyle=$
$\displaystyle{\mathcal{N}}(\mbox{\boldmath$g$}_{k};\;\mathcal{G}(\mbox{\boldmath$s$}_{n}),{\bf\Gamma}_{n}),$
(7)

with:

${\mathcal{N}}(\mbox{\boldmath$f$}_{m};\;\mathcal{F}(\mbox{\boldmath$s$}_{n}),{\bf\Sigma}_{n})=\frac{1}{(2\pi)^{r/2}|{\bf\Sigma}_{n}|^{1/2}}\exp\left(-\frac{1}{2}\|\mbox{\boldmath$f$}_{m}-\mathcal{F}(\mbox{\boldmath$s$}_{n})\|^{2}_{{\bf\Sigma}_{n}}\right),$
(8)

where the notation
$\|\mbox{\boldmath$v$}-\mbox{\boldmath$w$}\|^{2}_{\bf\Sigma}$ stands for the
Mahalanobis distance
$(\mbox{\boldmath$v$}-\mbox{\boldmath$w$})^{\top}{\bf\Sigma}^{-1}(\mbox{\boldmath$v$}-\mbox{\boldmath$w$})$
and ⊤ stands for the transpose of a matrix. The likelihoods of outliers are
taken as two uniform distributions:

$\displaystyle
P_{{N+1}}^{\mathbb{F}}(\mbox{\boldmath$f$}_{m})=\mathcal{U}(\mbox{\boldmath$f$}_{m};V),$
(9) $\displaystyle
P_{{N+1}}^{\mathbb{G}}(\mbox{\boldmath$g$}_{k})=\mathcal{U}(\mbox{\boldmath$g$}_{k};U),$
(10)

where $V$ and $U$ denote the respective support volumes. We also define the
prior probabilities $\pi=(\pi_{1},\ldots,\pi_{n},\ldots,\pi_{N+1})$ and
$\lambda=(\lambda_{1},\ldots,\lambda_{n},\ldots,\lambda_{N+1})$:

$\displaystyle\pi_{n}$ $\displaystyle=$ $\displaystyle P(A_{m}=n),\quad\forall
m=1\ldots M,$ (11) $\displaystyle\lambda_{n}$ $\displaystyle=$ $\displaystyle
P(B_{k}=n),\quad\forall k=1\ldots K.$ (12)

Therefore, $\mbox{\boldmath$f$}_{m}$ and $\mbox{\boldmath$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:

$\displaystyle P(\mbox{\boldmath$f$}_{m})$ $\displaystyle=$
$\displaystyle\sum\limits_{n=1}^{N}\pi_{n}{\mathcal{N}}(\mbox{\boldmath$f$}_{m};\;\mathcal{F}(\mbox{\boldmath$s$}_{n}),{\bf\Sigma}_{n})+\pi_{N+1}\mathcal{U}(\mbox{\boldmath$f$}_{m};V),$
(13) $\displaystyle P(\mbox{\boldmath$g$}_{k})$ $\displaystyle=$
$\displaystyle\sum\limits_{n=1}^{N}\lambda_{n}{\mathcal{N}}(\mbox{\boldmath$g$}_{k};\;\mathcal{G}(\mbox{\boldmath$s$}_{n}),{\bf\Gamma}_{n})+\lambda_{N+1}\mathcal{U}(\mbox{\boldmath$g$}_{k};U).$
(14)

The log-likelihood of the observed data can then be written as:

$\displaystyle{\mathcal{L}}(\mbox{\bf f},\mbox{\bf
g},\mbox{\boldmath$\theta$})$ $\displaystyle=$
$\displaystyle\sum\limits_{m=1}^{M}\log\left(\sum\limits_{n=1}^{N}\pi_{n}{\mathcal{N}}(\mbox{\boldmath$f$}_{m};\;\mathcal{F}(\mbox{\boldmath$s$}_{n}),{\bf\Sigma}_{n})+\pi_{N+1}\mathcal{U}(\mbox{\boldmath$f$}_{m};V)\right)+$
(15) $\displaystyle+$
$\displaystyle\sum\limits_{k=1}^{K}\log\left(\sum\limits_{n=1}^{N}\lambda_{n}{\mathcal{N}}(\mbox{\boldmath$g$}_{k};\;\mathcal{G}(\mbox{\boldmath$s$}_{n}),{\bf\Gamma}_{n})+\lambda_{N+1}\mathcal{U}(\mbox{\boldmath$g$}_{k};U)\right)$

where:

$\mbox{\boldmath$\theta$}=\\{\pi_{1},\ldots,\pi_{N},\pi_{N+1},\lambda_{1},\ldots,\lambda_{N},\lambda_{N+1},\mbox{\boldmath$s$}_{1},\ldots,\mbox{\boldmath$s$}_{N},{\bf\Sigma}_{1},\ldots,{\bf\Sigma}_{N},{\bf\Gamma}_{1},\ldots,{\bf\Gamma}_{N}\\}$
(16)

denotes the set of all unknown parameters to be estimated using a maximum
likelihood principle.

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.

The graphical representation of our conjugate mixture model is shown in Figure
1. We adopted the graphical notation introduced in (Bishop, 2006) 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 $\mbox{\boldmath$s$}_{1},\ldots\mbox{\boldmath$s$}_{N}$ shown in
between the two plates.

##  3 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 log-likelihood (15) is difficult to achieve. The expectation-
maximization (EM) algorithm (Dempster et al., 1977; McLachlan and Krishnan,
1996) 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 $\theta$ and $\tilde{\mbox{\boldmath$\theta$}}$:

$\displaystyle{\mathcal{L}}(\mbox{\bf f},\mbox{\bf
g},\mbox{\boldmath$\theta$})$ $\displaystyle=$ $\displaystyle
Q(\mbox{\boldmath$\theta$},\tilde{\mbox{\boldmath$\theta$}})+H(\mbox{\boldmath$\theta$},\tilde{\mbox{\boldmath$\theta$}}),$
(17) $\displaystyle\mbox{with}\quad
Q(\mbox{\boldmath$\theta$},\tilde{\mbox{\boldmath$\theta$}})$ $\displaystyle=$
$\displaystyle\mathbb{E}[\log P(\mbox{\bf f},\mbox{\bf
g},\mbox{\boldmath$A$},\mbox{\boldmath$B$};\mbox{\boldmath$\theta$})\;|\;\mbox{\bf
f},\mbox{\bf g};\tilde{\mbox{\boldmath$\theta$}}],$ (18)
$\displaystyle\mbox{and}\quad
H(\mbox{\boldmath$\theta$},\tilde{\mbox{\boldmath$\theta$}})$ $\displaystyle=$
$\displaystyle-\mathbb{E}[\log
P(\mbox{\boldmath$A$},\mbox{\boldmath$B$}\;|\;\mbox{\bf f},\mbox{\bf
g};\mbox{\boldmath$\theta$})|\mbox{\bf f},\mbox{\bf
g};\tilde{\mbox{\boldmath$\theta$}}],$ (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 $\mbox{\boldmath$\theta$}^{{(q)}}$ of the
parameters, compute the conditional expectation with respect to variables $A$
and $B$:

$Q(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{{(q)}})=\sum\limits_{\mbox{\boldmath$a$}\in\\{1\ldots{N+1}\\}^{M}}\sum\limits_{\mbox{\boldmath$b$}\in\\{1\ldots{N+1}\\}^{K}}P(\mbox{\boldmath$a$},\mbox{\boldmath$b$}|\mbox{\bf
f},\mbox{\bf g};\;\mbox{\boldmath$\theta$}^{{(q)}})\;\log P(\mbox{\bf
f},\mbox{\bf
g},\mbox{\boldmath$a$},\mbox{\boldmath$b$};\;\mbox{\boldmath$\theta$})$ (20)

  * •

Maximization. Update the parameter set $\mbox{\boldmath$\theta$}^{{(q)}}$ by
maximizing (20) with respect to $\theta$:

$\mbox{\boldmath$\theta$}^{{(q+1)}}=\operatornamewithlimits{argmax}\limits_{\mbox{\boldmath$\theta$}}Q(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{{(q)}})$
(21)

It is well known that the EM algorithm increases the target function
${\mathcal{L}}(\mbox{\bf f},\mbox{\bf g},\mbox{\boldmath$\theta$})$ in (15),
i.e., the sequence of estimates
$\\{\mbox{\boldmath$\theta$}^{{(q)}}\\}_{q\in\mathbb{N}}$ satisfies
${\mathcal{L}}(\mbox{\bf f},\mbox{\bf
g},\mbox{\boldmath$\theta$}^{{(q+1)}})\geq{\mathcal{L}}(\mbox{\bf f},\mbox{\bf
g},\mbox{\boldmath$\theta$}^{{(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 yields Generalized EM
(GEM) procedures (McLachlan and Krishnan, 1996) (see (Boyles, 1983) for a
result on the convergence of this class of algorithms). The GEM algorithm
searches for some $\mbox{\boldmath$\theta$}^{{(q+1)}}$ such that
$Q(\mbox{\boldmath$\theta$}^{{(q+1)}},\mbox{\boldmath$\theta$}^{{(q)}})\geq
Q(\mbox{\boldmath$\theta$}^{{(q)}},\mbox{\boldmath$\theta$}^{{(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.

###  3.1 The Expectation Step

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

$Q(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{{(q)}})=Q_{\mathcal{F}}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{{(q)}})+Q_{\mathcal{G}}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{{(q)}}),$
(22)

with

$\displaystyle
Q_{\mathcal{F}}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{{(q)}})$
$\displaystyle=\sum\limits_{m=1}^{M}\sum\limits_{n=1}^{N+1}\alpha_{mn}^{{(q)}}\log\big{(}\pi_{n}P(\mbox{\boldmath$f$}_{m}|A_{m}=n;\;\mbox{\boldmath$\theta$})\big{)},$
(23) $\displaystyle\quad
Q_{\mathcal{G}}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{{(q)}})$
$\displaystyle=\sum\limits_{k=1}^{K}\sum\limits_{n=1}^{N+1}\beta_{kn}^{{(q)}}\log\big{(}\lambda_{n}P(\mbox{\boldmath$g$}_{k}|B_{k}=n;\;\mbox{\boldmath$\theta$})\big{)},$
(24)

where $\alpha_{mn}^{{(q)}}$ and $\beta_{kn}^{{(q)}}$ denote the posterior
probabilities
$\alpha_{mn}^{{(q)}}=P(A_{m}=n|\mbox{\boldmath$f$}_{m};\mbox{\boldmath$\theta$}^{{(q)}})$
and
$\beta_{kn}^{{(q)}}=P(B_{k}=n|\mbox{\boldmath$g$}_{k};\mbox{\boldmath$\theta$}^{{(q)}})$.
Their expressions can be derived straightforwardly from Bayes’ theorem,
$\forall n=1\ldots N$:

$\displaystyle\alpha_{mn}^{{(q)}}=$
$\displaystyle\frac{\pi_{n}^{{(q)}}{\mathcal{N}}(\mbox{\boldmath$f$}_{m};\mathcal{F}(\mbox{\boldmath$s$}_{n}^{{(q)}}),{\bf\Sigma}_{n}^{{(q)}})}{\sum\limits_{i=1}^{N}\pi_{i}^{{(q)}}{\mathcal{N}}(\mbox{\boldmath$f$}_{m};\mathcal{F}(\mbox{\boldmath$s$}_{i}^{{(q)}}),{\bf\Sigma}_{i}^{{(q)}})+V^{-1}\pi_{N+1}^{{(q)}}},$
(25) $\displaystyle\quad\beta_{kn}^{{(q)}}=$
$\displaystyle\frac{\lambda_{n}^{{(q)}}{\mathcal{N}}(\mbox{\boldmath$g$}_{k};\;\mathcal{G}(\mbox{\boldmath$s$}_{n}^{{(q)}}),{\bf\Gamma}_{n}^{{(q)}})}{\sum\limits_{i=1}^{N}\lambda_{i}^{{(q)}}{\mathcal{N}}(\mbox{\boldmath$g$}_{k};\mathcal{G}(\mbox{\boldmath$s$}_{i}^{{(q)}}),{\bf\Gamma}_{i}^{{(q)}})+U^{-1}\lambda^{(q)}_{N+1}}.$
(26)

and $\alpha_{m,{N+1}}^{{(q)}}=1-\sum\limits_{n=1}^{N}\alpha_{mn}^{{(q)}}$ and
$\beta_{k,{N+1}}^{{(q)}}=1-\sum\limits_{n=1}^{N}\beta_{kn}^{{(q)}}$. Using
(6)-(10) the expressions above further lead to:

$\displaystyle
Q_{\mathcal{F}}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{{(q)}})=$
$\displaystyle-\frac{1}{2}\sum\limits_{m=1}^{M}\sum\limits_{n=1}^{N}\alpha_{mn}^{{(q)}}\left(\|\mbox{\boldmath$f$}_{m}-\mathcal{F}(\mbox{\boldmath$s$}_{n})\|^{2}_{{\bf\Sigma}_{n}}+\log((2\pi)^{r}|{\bf\Sigma}_{n}|\pi_{n}^{-2})\right)-$
$\displaystyle-\frac{1}{2}\sum\limits_{m=1}^{M}\alpha_{m,{N+1}}^{{(q)}}\log(V^{2}\pi_{N+1}^{-2}),$
(27) $\displaystyle
Q_{\mathcal{G}}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{{(q)}})=$
$\displaystyle-\frac{1}{2}\sum\limits_{k=1}^{K}\sum\limits_{n=1}^{N}\beta_{kn}^{{(q)}}\left(\|\mbox{\boldmath$g$}_{k}-\mathcal{G}(\mbox{\boldmath$s$}_{n})\|^{2}_{{\bf\Gamma}_{n}}+\log((2\pi)^{p}|{\bf\Gamma}_{n}|\lambda_{n}^{-2})\right)-$
$\displaystyle-\frac{1}{2}\sum\limits_{k=1}^{K}\beta_{k,{N+1}}^{{(q)}}\log(U^{2}\lambda_{N+1}^{-2}).$
(28)

###  3.2 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
$\forall n=1,\ldots,{N+1}$:

$\displaystyle\pi_{n}^{{(q+1)}}$ $\displaystyle=$
$\displaystyle\frac{1}{M}\sum\limits_{m=1}^{M}\alpha_{mn}^{{(q)}},$ (29)
$\displaystyle\lambda_{n}^{{(q+1)}}$ $\displaystyle=$
$\displaystyle\frac{1}{K}\sum\limits_{k=1}^{K}\beta_{kn}^{{(q)}}.$ (30)

The covariance matrices are governed by the tying parameters
$\mbox{\boldmath$s$}_{n}^{{(q+1)}}\in\mathbb{S}$ through the functions
$\mathcal{F}$ and $\mathcal{G}$, $\forall n=1,\ldots,N$:

$\displaystyle{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$}_{n}^{{(q+1)}})$
$\displaystyle=$
$\displaystyle\frac{1}{\sum\limits_{m=1}^{M}\alpha_{mn}^{{(q)}}}\sum\limits_{m=1}^{M}\alpha_{mn}^{{(q)}}(\mbox{\boldmath$f$}_{m}-\mathcal{F}(\mbox{\boldmath$s$}_{n}^{{(q+1)}}))(\mbox{\boldmath$f$}_{m}-\mathcal{F}(\mbox{\boldmath$s$}_{n}^{{(q+1)}}))^{\top},$
(31)
$\displaystyle{\bf\Gamma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$}_{n}^{{(q+1)}})$
$\displaystyle=$
$\displaystyle\frac{1}{\sum\limits_{k=1}^{K}\beta_{kn}^{{(q)}}}\sum\limits_{k=1}^{K}\beta_{kn}^{{(q)}}(\mbox{\boldmath$g$}_{k}-\mathcal{G}(\mbox{\boldmath$s$}_{n}^{{(q+1)}}))(\mbox{\boldmath$g$}_{k}-\mathcal{G}(\mbox{\boldmath$s$}_{n}^{{(q+1)}}))^{\top}.$
(32)

For every $n=1,\ldots,N$, $\mbox{\boldmath$s$}_{n}^{{(q+1)}}$ is the parameter
vector such that:

$\mbox{\boldmath$s$}_{n}^{{(q+1)}}=\operatornamewithlimits{argmax}\limits_{\mbox{\boldmath$s$}}Q_{n}^{{(q)}}(\mbox{\boldmath$s$}),$
(33)

where

$\displaystyle Q_{n}^{{(q)}}(\mbox{\boldmath$s$})=$ $\displaystyle-$
$\displaystyle\sum\limits_{m=1}^{M}\alpha_{mn}^{{(q)}}(\|\mbox{\boldmath$f$}_{m}-\mathcal{F}(\mbox{\boldmath$s$})\|^{2}_{{\bf\Sigma}_{n}(\mbox{\boldmath$s$})}+\log|{\bf\Sigma}_{n}(\mbox{\boldmath$s$})|)-$
(34) $\displaystyle-$
$\displaystyle\sum\limits_{k=1}^{K}\beta_{kn}^{{(q)}}(\|\mbox{\boldmath$g$}_{k}-\mathcal{G}(\mbox{\boldmath$s$})\|^{2}_{{\bf\Gamma}_{n}(\mbox{\boldmath$s$})}+\log|{\bf\Gamma}_{n}(\mbox{\boldmath$s$})|).$

We stress that the covariances ${\bf\Sigma}_{n}(\mbox{\boldmath$s$})$ and
${\bf\Gamma}_{n}(\mbox{\boldmath$s$})$ in (31) and (32) are considered as
functions of $\mbox{\boldmath$s$}\in\mathbb{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, $\mathcal{F}$ and $\mathcal{G}$ are non-linear
transformations and hence there is no simple closed-form expression for the
estimation of the tying parameters.

###  3.3 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. 1.

Apply procedure Initialize to initialize the parameter vector:
$\mbox{\boldmath$\theta$}^{(0)}=\\{\pi_{1}^{(0)},\ldots,\pi_{N+1}^{(0)},\lambda_{1}^{(0)},\ldots,\lambda_{N+1}^{(0)},\mbox{\boldmath$s$}_{1}^{(0)},\ldots,\mbox{\boldmath$s$}_{N}^{(0)},{\bf\Sigma}_{1}^{(0)},\ldots,{\bf\Sigma}_{N}^{(0)},{\bf\Gamma}_{1}^{(0)},\ldots,{\bf\Gamma}_{N}^{(0)}\\}$;

  2. 2.

E step: compute $Q(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{{(q)}})$
using equations (25) to (28);

  3. 3.

M step: estimate $\mbox{\boldmath$\theta$}^{{(q+1)}}$ using the following sub-
steps:

    1. (a)

The priors. Compute $\pi_{1}^{{(q+1)}},\ldots,\pi_{N+1}^{{(q+1)}}$ and
$\lambda_{1}^{{(q+1)}},\ldots,\lambda_{N+1}^{{(q+1)}}$ using (29) and (30);

    2. (b)

The tying parameters. For each $n=1\ldots N$:

       * •

Apply procedure Choose to determine an initial value, denoted by
$\tilde{\mbox{\boldmath$s$}}_{n}^{(0)}$, as proposed in Section 5;

       * •

Apply procedure Local Search to each $Q_{n}^{{(q)}}(\mbox{\boldmath$s$})$ as
defined in (34) starting from $\tilde{\mbox{\boldmath$s$}}_{n}^{(0)}$ and set
the result to $\mbox{\boldmath$s$}_{n}^{{(q+1)}}$ using the eq. (35) specified
below;

    3. (c)

The covariance matrices. For every $n=1\ldots N$, use (31) and (32) to compute
${\bf\Sigma}_{n}^{{(q+1)}}$ and ${\bf\Gamma}_{n}^{{(q+1)}}$;

  4. 4.

Check for convergence: Terminate, otherwise go to Step 2;

  5. 5.

Apply procedure Select, use (62) specified below to determine the best $N$;

This algorithm uses the following procedures:

  * •

Initialize: this procedure aims at providing the initial parameter values
$\mbox{\boldmath$\theta$}^{(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 single-space 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
$\tilde{\mbox{\boldmath$s$}}_{1}^{(0)},\ldots,\tilde{\mbox{\boldmath$s$}}_{N}^{(0)}$
which are likely to be close to the global maxima of the functions
$Q_{n}^{{(q)}}(\mbox{\boldmath$s$})$ in (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
(Zhigljavsky and Žilinskas, 2008) and propose different variants in Section 5.

  * •

Local Search: an important requirement of this procedure is that it finds a
local maximum of the $Q_{n}^{{(q)}}(\mbox{\boldmath$s$})$’s starting from any
arbitrary point in $\mathbb{S}$. In this work, we will consider procedures
that consist in iterating a local update of the form ($\nu$ is the iteration
index):

$\tilde{\mbox{\boldmath$s$}_{n}}^{(\nu+1)}=\tilde{\mbox{\boldmath$s$}_{n}}^{(\nu)}+\mbox{\bf
H}_{n}^{(q,\nu)}\mbox{\boldmath$\nabla$}Q_{n}^{{(q)}}(\tilde{\mbox{\boldmath$s$}_{n}}^{(\nu)}),$
(35)

with $\mbox{\bf H}_{n}^{(q,\nu)}$ being a positive definite matrix that may
vary with $\nu$. When the gradient
$\mbox{\boldmath$\nabla$}Q_{n}^{{(q)}}(\mbox{\boldmath$s$})$ is Lipschitz
continuous with some constant $L_{n}^{{(q)}}$, an appropriate choice that
guarantees the increase of
$Q_{n}^{{(q)}}(\tilde{\mbox{\boldmath$s$}}^{(\nu)})$ at each iteration $\nu$,
is to choose $\mbox{\bf H}_{n}^{(q,\nu)}$ such that it verifies $\|\mbox{\bf
H}_{n}^{\,(q,\nu)}\|\leq 2/L_{n}^{{(q)}}$.

Different choices for $\mbox{\bf H}_{n}^{(q,\nu)}$ are possible and they
correspond to different optimization methods that belong, in general, to the
variable metric class. For example $\mbox{\bf
H}_{n}^{(q,\nu)}=\frac{2}{L_{n}^{{(q)}}}\mbox{\bf I}$ leads to gradient
ascent, while taking $\mbox{\bf H}_{n}^{(q,\nu)}$ 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.

##  4 Analysis of the Local Search Procedure

Each instance of (34) for $n=1,\ldots,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
$Q_{n}^{{(q)}}(\mbox{\boldmath$s$})$ in (34) is increased. We start by
rewriting (34) more conveniently in order to perform the optimization with
respect to $\mbox{\boldmath$s$}\in\mathbb{S}$. To simplify the notation, the
iteration index $q$ is sometimes omitted. We simply write
$Q_{n}(\mbox{\boldmath$s$})$ for $Q_{n}^{{(q)}}(\mbox{\boldmath$s$})$.

Let $\bar{\alpha}_{n}=\sum_{m=1}^{M}\alpha_{mn}^{{(q)}}$ and
$\bar{\beta}_{n}=\sum_{k=1}^{K}\beta_{kn}^{{(q)}}$ denote the average object
weights in each one of the two modalities. We introduce
$\alpha_{n}=\bar{\alpha}_{n}^{-1}(\alpha^{(q)}_{1n},\ldots,\alpha^{(q)}_{Mn})$
and $\beta_{n}=\bar{\beta}_{n}^{-1}(\beta^{(q)}_{1n},\ldots,\beta^{(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
$\\{\mbox{\boldmath$f$}_{1},\ldots,\mbox{\boldmath$f$}_{m},\ldots,\mbox{\boldmath$f$}_{M}\\}$
and
$\\{\mbox{\boldmath$g$}_{1},\ldots,\mbox{\boldmath$g$}_{k},\ldots,\mbox{\boldmath$g$}_{K}\\}$.
It follows that the expressions for the optimal variances (31) and (32) as
functions of $s$, can be rewritten as:

$\displaystyle{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})=$
$\displaystyle\mathbb{E}_{\alpha_{n}}[\big{(}\mbox{\boldmath$F$}-\mathcal{F}(\mbox{\boldmath$s$})\big{)}\big{(}\mbox{\boldmath$F$}-\mathcal{F}(\mbox{\boldmath$s$})\big{)}^{\top}],$
(36) $\displaystyle{\bf\Gamma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})=$
$\displaystyle\mathbb{E}_{\beta_{n}}[\big{(}\mbox{\boldmath$G$}-\mathcal{G}(\mbox{\boldmath$s$})\big{)}\big{(}\mbox{\boldmath$G$}-\mathcal{G}(\mbox{\boldmath$s$})\big{)}^{\top}],$
(37)

where $\mathbb{E}_{\alpha_{n}}$ and $\mathbb{E}_{\beta_{n}}$ denote the
expectations with respect to the distributions $\alpha_{n}$ and $\beta_{n}$.
Using some standard projection formula, it follows that the covariances are:

$\displaystyle{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})=$
$\displaystyle\mbox{\bf
V}_{f}+\mbox{\boldmath$v$}_{f}\mbox{\boldmath$v$}_{f}^{\top},$ (38)
$\displaystyle{\bf\Gamma}^{{(q+1)}}(\mbox{\boldmath$s$})=$
$\displaystyle\mbox{\bf
V}_{g}+\mbox{\boldmath$v$}_{g}\mbox{\boldmath$v$}_{g}^{\top},$ (39)

where $\mbox{\bf V}_{f}$ and $\mbox{\bf V}_{g}$ are the covariance matrices of
$F$ and $G$ respectively under distributions $\alpha_{n}$ and $\beta_{n}$, and
$\mbox{\boldmath$v$}_{f}$ and $\mbox{\boldmath$v$}_{g}$ are vectors defined
by:

$\displaystyle\mbox{\boldmath$v$}_{f}=\mathbb{E}_{\alpha_{n}}[\mbox{\boldmath$F$}]-\mathcal{F}(\mbox{\boldmath$s$}),$
(40)
$\displaystyle\mbox{\boldmath$v$}_{g}=\mathbb{E}_{\beta_{n}}[\mbox{\boldmath$G$}]-\mathcal{G}(\mbox{\boldmath$s$}).$
(41)

For convenience we omit the index $n$ for $\mbox{\bf V}_{f}$, $\mbox{\bf
V}_{g}$, $\mbox{\boldmath$v$}_{f}$ and $\mbox{\boldmath$v$}_{g}$. Let
$\bar{\mbox{\boldmath$f$}}_{n}=\mathbb{E}_{\alpha_{n}}[\mbox{\boldmath$F$}]$
and
$\bar{\mbox{\boldmath$g$}}_{n}=\mathbb{E}_{\beta_{n}}[\mbox{\boldmath$G$}]$.
This yields:

$\displaystyle\bar{\mbox{\boldmath$f$}}_{n}$
$\displaystyle=\bar{\alpha}_{n}^{-1}\sum\limits_{m=1}^{M}\alpha_{mn}^{{(q)}}\mbox{\boldmath$f$}_{m},$
(42) $\displaystyle\bar{\mbox{\boldmath$g$}}_{n}$
$\displaystyle=\bar{\beta}_{n}^{-1}\sum\limits_{k=1}^{K}\beta_{kn}^{{(q)}}\mbox{\boldmath$g$}_{k},$
(43) $\displaystyle\mbox{\bf V}_{f}$
$\displaystyle=\bar{\alpha}_{n}^{-1}\sum\limits_{m=1}^{M}\alpha_{mn}^{{(q)}}\mbox{\boldmath$f$}_{m}\mbox{\boldmath$f$}_{m}^{\top}-\bar{\mbox{\boldmath$f$}}_{n}\bar{\mbox{\boldmath$f$}}_{n}^{\top},$
(44) $\displaystyle\mbox{\bf V}_{g}$
$\displaystyle=\bar{\beta}_{n}^{-1}\sum\limits_{k=1}^{K}\beta_{kn}^{{(q)}}\mbox{\boldmath$g$}_{k}\mbox{\boldmath$g$}_{k}^{\top}-\bar{\mbox{\boldmath$g$}}_{n}\bar{\mbox{\boldmath$g$}}_{n}^{\top}.$
(45)

Next we derive a simplified expression for $Q_{n}(\mbox{\boldmath$s$})$ in
(34) in order to investigate its properties. Notice that one can write (34) as
the sum
$Q_{n}(\mbox{\boldmath$s$})=Q_{n,\mathcal{F}}(\mbox{\boldmath$s$})+Q_{n,\mathcal{G}}(\mbox{\boldmath$s$})$,
with:

$Q_{n,\mathcal{F}}(\mbox{\boldmath$s$})=-\sum\limits_{m=1}^{M}\alpha_{mn}^{{(q)}}(\|\mbox{\boldmath$f$}_{m}-\mathcal{F}(\mbox{\boldmath$s$})\|^{2}_{{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})}+\log|{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})|),$
(46)

and a similar expression for $Q_{n,\mathcal{G}}(\mbox{\boldmath$s$})$. Eq.
(46) can be written:

$Q_{n,\mathcal{F}}(\mbox{\boldmath$s$})=-\bar{\alpha}_{n}(\mathbb{E}_{\alpha_{n}}[(\mbox{\boldmath$F$}-\mathcal{F}(\mbox{\boldmath$s$}))^{\top}{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})^{-1}(\mbox{\boldmath$F$}-\mathcal{F}(\mbox{\boldmath$s$}))]+\log|{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})|).$
(47)

The first term of (47) can be further divided into two terms:

$\displaystyle\mathbb{E}_{\alpha_{n}}[(\mbox{\boldmath$F$}-\mathcal{F}(\mbox{\boldmath$s$}))^{\top}{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})^{-1}(\mbox{\boldmath$F$}-\mathcal{F}(\mbox{\boldmath$s$}))]=$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\alpha_{n}}[(\mbox{\boldmath$F$}-\bar{\mbox{\boldmath$f$}}_{n})^{\top}{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})^{-1}(\mbox{\boldmath$F$}-\bar{\mbox{\boldmath$f$}}_{n})]+\mbox{\boldmath$v$}_{f}^{\top}{{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})}^{-1}\mbox{\boldmath$v$}_{f}.$
(48)

The Sherman-Morrison formula applied to (38) leads to

${\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})^{-1}=\mbox{\bf
V}_{f}^{-1}-\mbox{\bf
V}_{f}^{-1}\mbox{\boldmath$v$}_{f}\mbox{\boldmath$v$}_{f}^{\top}\mbox{\bf
V}_{f}^{-1}/(1+D_{n,\mathcal{F}}(\mbox{\boldmath$s$})),$ (49)

with:

$D_{n,\mathcal{F}}(\mbox{\boldmath$s$})=\|\mathcal{F}(\mbox{\boldmath$s$})-\bar{\mbox{\boldmath$f$}}_{n}\|_{\mbox{\bf
V}_{f}}^{2}.$ (50)

It follows that (48) can be written as the sum of:

$\mathbb{E}_{\alpha_{n}}[(\mbox{\boldmath$F$}-\bar{\mbox{\boldmath$f$}}_{n})^{\top}{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})^{-1}(\mbox{\boldmath$F$}-\bar{\mbox{\boldmath$f$}}_{n})]=C_{f}-\frac{D_{n,\mathcal{F}}(\mbox{\boldmath$s$})}{1+D_{n,\mathcal{F}}(\mbox{\boldmath$s$})},$
(51)

and of

$\mbox{\boldmath$v$}_{f}^{\top}{{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})}^{-1}\mbox{\boldmath$v$}_{f}=\frac{D_{n,\mathcal{F}}(\mbox{\boldmath$s$})}{1+D_{n,\mathcal{F}}(\mbox{\boldmath$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:

$\displaystyle\log|{\bf\Sigma}_{n}^{{(q+1)}}(\mbox{\boldmath$s$})|$
$\displaystyle=\log|\mbox{\bf
V}_{f}+\mbox{\boldmath$v$}_{f}\mbox{\boldmath$v$}_{f}^{\top}|=\log|\mbox{\bf
V}_{f}|+\log(1+\mbox{\boldmath$v$}_{f}^{\top}\mbox{\bf
V}_{f}^{-1}\mbox{\boldmath$v$}_{f})=$ $\displaystyle=\log|\mbox{\bf
V}_{f}|+\log(1+D_{n,\mathcal{F}}(\mbox{\boldmath$s$})).$ (53)

It follows that there is only one term depending on $s$ in (47):

$Q_{n,\mathcal{F}}(\mbox{\boldmath$s$})=-\bar{\alpha}_{n}\left(C_{f}+\log|\mbox{\bf
V}_{f}|+\log(1+D_{n,\mathcal{F}}(\mbox{\boldmath$s$}))\right).$ (54)

Repeating the same derivation for the second sensorial modality we obtain the
following equivalent form of (34):

$Q_{n}(\mbox{\boldmath$s$})=-\bar{\alpha}_{n}\log(1+D_{n,\mathcal{F}}(\mbox{\boldmath$s$}))-\bar{\beta}_{n}\log(1+D_{n,\mathcal{G}}(\mbox{\boldmath$s$}))+C,$
(55)

where $C$ is some constant not depending on $s$.

Using this form of $Q_{n}(\mbox{\boldmath$s$})$, we can now investigate the
properties of its gradient
$\mbox{\boldmath$\nabla$}Q_{n}(\mbox{\boldmath$s$})$. It appears that under
some regularity assumptions on $\mathcal{F}$ and $\mathcal{G}$, the gradient
$\mbox{\boldmath$\nabla$}Q_{n}(\mbox{\boldmath$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
$\|\mbox{\bf V}\|=\sup\limits_{\|\mbox{\boldmath$v$}\|=1}\|\mbox{\bf
V}\mbox{\boldmath$v$}\|$. For simplicity, we further omit the index $n$.

######  Lemma 1.

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

$\varphi(\mbox{\boldmath$v$})=\|\mbox{\bf
V}\mbox{\boldmath$v$}\|/(1+\mbox{\boldmath$v$}^{\top}\mbox{\bf
V}\mbox{\boldmath$v$})$

is bounded by $\varphi(\mbox{\boldmath$v$})\leq C_{\varphi}(\mbox{\bf V})$
with $C_{\varphi}(\mbox{\bf V})=\sqrt{\|\mbox{\bf V}\|}/2$ and is Lipschitz
continuous:

$\forall\mbox{\boldmath$v$},\tilde{\mbox{\boldmath$v$}}\quad\|\varphi(\mbox{\boldmath$v$})-\varphi(\tilde{\mbox{\boldmath$v$}})\|\leq
L_{\varphi}(\mbox{\bf V})\|\mbox{\boldmath$v$}-\tilde{\mbox{\boldmath$v$}}\|,$

where $L_{\varphi}(\mbox{\bf V})=\|\mbox{\bf V}\|(1+\mu(\mbox{\bf V})/2)$ is
the Lipschitz constant and $\mu(\mbox{\bf V})=\|\mbox{\bf V}\|\|\mbox{\bf
V}^{-1}\|$ is the condition number of V.

Proof:  We start by introducing $\mbox{\boldmath$w$}=\mbox{\bf
V}\mbox{\boldmath$v$}$ so that
$\varphi(\mbox{\boldmath$v$})=\tilde{\varphi}(\mbox{\boldmath$w$})=\|\mbox{\boldmath$w$}\|/(1+\mbox{\boldmath$w$}^{\top}\mbox{\bf
V}^{-1}\mbox{\boldmath$w$}).$ As soon as $\mbox{\boldmath$w$}^{\top}\mbox{\bf
V}^{-1}\mbox{\boldmath$w$}\geq\lambda_{\mathrm{min}}\|\mbox{\boldmath$w$}\|^{2}$
(where we denoted by $\lambda_{\mathrm{min}}$ the smallest eigenvalue of
$\mbox{\bf V}^{-1}$, so that in fact $\lambda_{\mathrm{min}}=\|\mbox{\bf
V}\|^{-1}$), to find the maximum of $\tilde{\varphi}(\mbox{\boldmath$w$})$ we
should maximize the expression $t/(1+\lambda_{\mathrm{min}}t^{2})$ for
$t=\|\mbox{\boldmath$w$}\|\geq 0$. It is reached at the point
$t^{*}=\lambda_{\mathrm{min}}^{-1/2}$. Substituting this value into the
original expressions gives $\varphi(\mbox{\boldmath$v$})\leq\sqrt{\|\mbox{\bf
V}\|}/2$.

To compute the Lipschitz constant $L_{\varphi}$ we consider the derivative:

$\displaystyle\|\mbox{\boldmath$\nabla$}\tilde{\varphi}^{\prime}(\mbox{\boldmath$w$})\|$
$\displaystyle=\frac{\left\|(1+\mbox{\boldmath$w$}^{\top}\mbox{\bf
V}^{-1}\mbox{\boldmath$w$})\mbox{\boldmath$w$}-2\|\mbox{\boldmath$w$}\|^{2}\mbox{\bf
V}^{-1}\mbox{\boldmath$w$}\right\|}{\|\mbox{\boldmath$w$}\|(1+\mbox{\boldmath$w$}^{\top}\mbox{\bf
V}^{-1}\mbox{\boldmath$w$})^{2}}\leq 1+\frac{2\|\mbox{\bf
V}^{-1}\|\|\mbox{\boldmath$w$}\|^{2}}{(1+\mbox{\boldmath$w$}^{\top}\mbox{\bf
V}^{-1}\mbox{\boldmath$w$})^{2}},$

from where we find that
$\|\mbox{\boldmath$\nabla$}\tilde{\varphi}^{\prime}(\mbox{\boldmath$w$})\|\leq
1+\mu(\mbox{\bf V})/2$, and so $L_{\varphi}=\|\mbox{\bf V}\|(1+\mu(\mbox{\bf
V})/2).\quad\blacksquare$

This lemma yields the following main result for the gradient
$\mbox{\boldmath$\nabla$}Q$:

######  Theorem 1.

Assume functions $\mathcal{F}$ and $\mathcal{G}$ and their derivatives
$\mathcal{F}^{\prime}$ and $\mathcal{G}^{\prime}$ are Lipschitz continuous
with constants $L_{\mathcal{F}}$, $L_{\mathcal{G}}$,
$L^{\prime}_{\mathcal{F}}$ and $L^{\prime}_{\mathcal{G}}$ respectively. Then
the gradient $\mbox{\boldmath$\nabla$}Q$ is bounded and Lipschitz continuous
with some constant $L$.

Proof:  From (55) the gradient $\mbox{\boldmath$\nabla$}Q$ can be written as:

$\displaystyle\mbox{\boldmath$\nabla$}Q(\mbox{\boldmath$s$})$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\nabla$}Q_{\mathcal{F}}(\mbox{\boldmath$s$})+\mbox{\boldmath$\nabla$}Q_{\mathcal{G}}(\mbox{\boldmath$s$})=$
(56) $\displaystyle=$
$\displaystyle\frac{2\bar{\alpha}{\mathcal{F}^{\prime}}^{\top}(\mbox{\boldmath$s$})\mbox{\bf
V}_{f}^{-1}(\bar{\mbox{\boldmath$f$}}-\mathcal{F}(\mbox{\boldmath$s$}))}{1+D_{\mathcal{F}}(\mbox{\boldmath$s$})}+\frac{2\bar{\beta}{\mathcal{G}^{\prime}}^{\top}(\mbox{\boldmath$s$})\mbox{\bf
V}_{g}^{-1}(\bar{\mbox{\boldmath$g$}}-\mathcal{G}(\mbox{\boldmath$s$}))}{1+D_{\mathcal{G}}(\mbox{\boldmath$s$})}.$

It follows from Lemma 1 that
$\left\|\mbox{\boldmath$\nabla$}Q_{\mathcal{F}}(\mbox{\boldmath$s$})\right\|\leq
2L_{\mathcal{F}}\bar{\alpha}C_{\varphi}(\mbox{\bf V}_{f}^{-1})\;$ and
$\;\left\|\mbox{\boldmath$\nabla$}Q_{\mathcal{G}}(\mbox{\boldmath$s$})\right\|\leq
2L_{\mathcal{G}}\bar{\beta}C_{\varphi}(\mbox{\bf V}_{g}^{-1})$ . The norm of
the gradient is then bounded by:

$\left\|\mbox{\boldmath$\nabla$}Q(\mbox{\boldmath$s$})\right\|\leq
2L_{\mathcal{F}}\bar{\alpha}C_{\varphi}(\mbox{\bf
V}_{f}^{-1})+2L_{\mathcal{G}}\bar{\beta}C_{\varphi}(\mbox{\bf V}_{g}^{-1}).$
(57)

Considering the norm
$\left\|\mbox{\boldmath$\nabla$}Q_{\mathcal{F}}(\mbox{\boldmath$s$})-\mbox{\boldmath$\nabla$}Q_{\mathcal{F}}(\tilde{\mbox{\boldmath$s$}})\right\|$,
we introduce
$\mbox{\boldmath$v$}_{1}=\bar{\mbox{\boldmath$f$}}-\mathcal{F}(\mbox{\boldmath$s$})$
and
$\mbox{\boldmath$v$}_{2}=\bar{\mbox{\boldmath$f$}}-\mathcal{F}(\tilde{\mbox{\boldmath$s$}})$.
Then we have:

$\displaystyle\left\|\mbox{\boldmath$\nabla$}Q_{\mathcal{F}}(\mbox{\boldmath$s$})-\mbox{\boldmath$\nabla$}Q_{\mathcal{F}}(\tilde{\mbox{\boldmath$s$}})\right\|$
$\displaystyle\leq
2\bar{\alpha}\left(\left\|\frac{(\mathcal{F}^{\prime}(\mbox{\boldmath$s$})-\mathcal{F}^{\prime}(\tilde{\mbox{\boldmath$s$}}))^{\top}\mbox{\bf
V}_{f}^{-1}\mbox{\boldmath$v$}_{1}}{1+\|\mbox{\boldmath$v$}_{1}\|^{2}_{\mbox{\bf
V}_{f}}}\right\|\right.+$
$\displaystyle+\left.\left\|\frac{{\mathcal{F}^{\prime}}^{\top}(\tilde{\mbox{\boldmath$s$}})\mbox{\bf
V}_{f}^{-1}\mbox{\boldmath$v$}_{2}}{1+\|\mbox{\boldmath$v$}_{2}\|^{2}_{\mbox{\bf
V}_{f}}}-\frac{{\mathcal{F}^{\prime}}^{\top}(\tilde{\mbox{\boldmath$s$}})\mbox{\bf
V}_{f}^{-1}\mbox{\boldmath$v$}_{1}}{1+\|\mbox{\boldmath$v$}_{1}\|^{2}_{\mbox{\bf
V}_{f}}}\right\|\right).$ (58)

Using Lemma 1 with $\mbox{\bf V}_{f}^{-1}$ we have:

$\left\|\mbox{\boldmath$\nabla$}Q_{\mathcal{F}}(\mbox{\boldmath$s$})-\mbox{\boldmath$\nabla$}Q_{\mathcal{F}}(\tilde{\mbox{\boldmath$s$}})\right\|\leq
2\bar{\alpha}\big{(}L^{\prime}_{\mathcal{F}}C_{\varphi}(\mbox{\bf
V}_{f}^{-1})+L^{2}_{\mathcal{F}}L_{\varphi}(\mbox{\bf
V}_{f}^{-1})\big{)}\|\mbox{\boldmath$s$}-\tilde{\mbox{\boldmath$s$}}\|.$

The same derivations can be performed for
$\mbox{\boldmath$\nabla$}Q_{\mathcal{G}}(\mbox{\boldmath$s$})$, so that
finally we get:

$\left\|\mbox{\boldmath$\nabla$}Q_{\mathcal{G}}(\mbox{\boldmath$s$})-\mbox{\boldmath$\nabla$}Q_{\mathcal{G}}(\tilde{\mbox{\boldmath$s$}})\right\|\leq
L\|\mbox{\boldmath$s$}-\tilde{\mbox{\boldmath$s$}}\|,$ (59)

where the Lipschitz constant is given by:

$L=2\bar{\alpha}\left(L^{\prime}_{\mathcal{F}}C_{\varphi}(\mbox{\bf
V}_{f}^{-1})+L^{2}_{\mathcal{F}}L_{\varphi}(\mbox{\bf
V}_{f}^{-1})\right)+2\bar{\beta}\left(L^{\prime}_{\mathcal{G}}C_{\varphi}(\mbox{\bf
V}_{g}^{-1})+L^{2}_{\mathcal{G}}L_{\varphi}(\mbox{\bf V}_{g}^{-1})\right).$
(60)

$\blacksquare$

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 ((Polyak, 1987)).

Let the function $Q:\mathbb{R}^{d}\to\mathbb{R}$ be differentiable on
$\mathbb{R}^{d}$ and its gradient $\mbox{\boldmath$\nabla$}Q$ be Lipschitz
continuous with constant $L$. Let the matrix H be positive definite, such that
$\|\mbox{\bf H}\|\leq\frac{2}{L}$. Then the sequence
$Q(\tilde{\mbox{\boldmath$s$}}^{(\nu)})$, defined by
$\tilde{\mbox{\boldmath$s$}}^{(\nu+1)}=\tilde{\mbox{\boldmath$s$}}^{(\nu)}+\mbox{\bf
H}\mbox{\boldmath$\nabla$}Q(\tilde{\mbox{\boldmath$s$}}^{(\nu)})$ is non-
decreasing.

This result shows that for any functions $\mathcal{F}$ and $\mathcal{G}$ that
verify the conditions of Theorem 1, using (35) with $\mbox{\bf
H}=\frac{2}{L}{\bf 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 $\nu$, if when considering the functions $\mathcal{F}$
and $\mathcal{G}$ on some restricted domain $\mathbb{S}^{(\nu)}$ there exist
smaller local Lipschitz constants $L^{(\nu)}_{\mathcal{F}}$,
$L^{(\nu)}_{\mathcal{G}}$, $L^{\prime(\nu)}_{\mathcal{F}}$ and
$L^{\prime(\nu)}_{\mathcal{G}}$, H can be set to $\mbox{\bf
H}=\frac{2}{L^{(\nu)}}{\bf I}$ with $L^{(\nu)}$ smaller than $L$. It follows
that
$\|\tilde{\mbox{\boldmath$s$}}^{(\nu+1)}-\tilde{\mbox{\boldmath$s$}}^{(\nu)}\|\leq\frac{2}{L^{(\nu)}}\|\mbox{\boldmath$\nabla$}Q(\tilde{\mbox{\boldmath$s$}}^{(\nu)})\|$,
which means that one can take the local constants, $L^{(\nu)}_{\mathcal{F}}$,
$L^{(\nu)}_{\mathcal{G}}$, $L^{\prime(\nu)}_{\mathcal{F}}$ and
$L^{\prime(\nu)}_{\mathcal{G}}$ if they are valid in the ball
$\mathbb{B}_{\rho^{(\nu)}}(\tilde{\mbox{\boldmath$s$}}^{(\nu)})$ with

$\rho^{(\nu)}=\frac{2}{L^{(\nu)}}\left(2L^{(\nu)}_{\mathcal{F}}\bar{\alpha}C_{\varphi}(\mbox{\bf
V}_{f}^{-1})+2L^{(\nu)}_{\mathcal{G}}\bar{\beta}C_{\varphi}(\mbox{\bf
V}_{g}^{-1})\right).$ (61)

##  5 Global Search and the Choose Procedure

Theorem 1 allows us to use the improved global random search techniques for
Lipschitz continuous functions (Zhigljavsky, 1991). 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 (Zhigljavsky, 1991). A proper choice of the initial value
$\tilde{\mbox{\boldmath$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 $\mathbb{S}$, according to some sequence of
distributions over $\mathbb{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 $\bar{\mbox{\boldmath$f$}}$ and $\bar{\mbox{\boldmath$g$}}$.
If there exists $\tilde{\mbox{\boldmath$s$}}^{(0)}$ such that
$\tilde{\mbox{\boldmath$s$}}^{(0)}\in\mathcal{F}^{-1}(\bar{\mbox{\boldmath$f$}})\cap\mathcal{G}^{-1}(\bar{\mbox{\boldmath$g$}})$,
then it is the global maximum and the M step solution is found. Otherwise, one
can sample $\mathbb{S}$ in the vicinity of the set
$\mathcal{F}^{-1}(\bar{\mbox{\boldmath$f$}})\cup\mathcal{G}^{-1}(\bar{\mbox{\boldmath$g$}})$
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 (Zhigljavsky, 1991). 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
$\tilde{\mbox{\boldmath$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.

##  6 Algorithm Initialization and the Initialize Procedure

In this section we focus on the problem of selecting the initial values
$\mbox{\boldmath$\theta$}^{(0)}$ for the model parameters. As it is often the
case with iterative optimization algorithms, the closer
$\mbox{\boldmath$\theta$}^{(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 initialization strategies, namely
the Observation Space Candidates (OSC) strategy and the Parameter Space
Candidates (PSC) strategy, that attempt to find a good initialization.

The Observation Space Candidates strategy consists in searching for cluster
centers in single modality spaces $\mathbb{F}$ and $\mathbb{G}$ to further map
them into the parameter space $\mathbb{S}$, and select the best candidates.
More specifically, we randomly select an observation $\mbox{\boldmath$f$}_{m}$
(or $\mbox{\boldmath$g$}_{k}$) and run the mean shift algorithm (Comaniciu and
Meer, 2002) in the corresponding space to find local modes of the
distribution, which are called candidates. The sets of candidate points
$\\{\hat{\mbox{\boldmath$f$}}_{i}\\}_{i\in I}$ and
$\\{\hat{\mbox{\boldmath$g$}}_{j}\\}_{j\in J}$ are further rarefied, that is
if
$\|\hat{\mbox{\boldmath$f$}}_{i_{1}}-\hat{\mbox{\boldmath$f$}}_{i_{2}}\|\leq\varepsilon_{\mathrm{f}}$
for some $i_{1}\neq i_{2}$ and for some threshold $\varepsilon>0$, we
eliminate one of these points. These rarefied sets are then mapped to
$\mathbb{S}$. If one of the observation space mappings, for example
$\mathcal{F}$, is non-injective, for each $\hat{\mbox{\boldmath$f$}}_{i}$ we
need to select a point
$\mbox{\boldmath$s$}_{i}\in\mathcal{F}^{-1}(\hat{\mbox{\boldmath$f$}}_{i})$
that is the best in some sense. We consider observations density in the other
observation spaces around an image of $\mbox{\boldmath$s$}_{i}$ as the
optimality measure of $\mbox{\boldmath$s$}_{i}$. This can be estimated through
calculation of the k-th nearest neighbour distance (k-NN) in the corresponding
observation space. The final step is to choose $N$ points out of these
candidates to initialize the cluster centers
$\\{\mbox{\boldmath$s$}_{1},\ldots,\mbox{\boldmath$s$}_{N}\\}$, so that the
inter-cluster distances are maximized. This can be done using, for example,
hierarchical clustering. The variances
${\bf\Sigma}_{1},\ldots,{\bf\Sigma}_{N}$ and
${\bf\Gamma}_{1},\ldots,{\bf\Gamma}_{N}$ are then calculated by standard
empirical variance formulas based on observations, that are closest to the
corresponding class center. The priors $\pi_{1},\ldots,\pi_{{N+1}}$ and
$\lambda_{1},\ldots,\lambda_{{N+1}}$ are set to be equal.

The Parameter Space Candidates strategy consists in mapping all the
observations to the parameter space $\mathbb{S}$, and performing subsequent
clustering in that space. More specifically, for every observation
$\mbox{\boldmath$f$}_{m}$ and $\mbox{\boldmath$g$}_{k}$ we find an optimal
point from the corresponding preimage
$\mathcal{F}^{-1}(\mbox{\boldmath$f$}_{m})$ and
$\mathcal{G}^{-1}(\mbox{\boldmath$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 $\mathbb{S}$ using the mean-shift algorithm, and initializing $N$
cluster centers $\\{\mbox{\boldmath$s$}_{1},\ldots,\mbox{\boldmath$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 10). 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.

##  7 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 (Schwarz,
1978). We consider the score function of the form

$\mathrm{BIC}_{N}=-2{\mathcal{L}}(\mbox{\bf f},\mbox{\bf
g},\hat{\mbox{\boldmath$\theta$}}_{N})+D_{N}\log(M+K),$ (62)

where $\hat{\mbox{\boldmath$\theta$}}_{N}$ is the ML estimate obtained by the
proposed EM algorithm, ${\mathcal{L}}(\mbox{\bf f},\mbox{\bf
g},\mbox{\boldmath$\theta$})$ is given by (15) and
$D_{N}=N\left(d+2+\frac{1}{2}(r^{2}+p^{2}+r+p)\right)$ 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(\mbox{\bf
f},\mbox{\bf g}|N=N_{0})$ because of the Hessian matrix of
${\mathcal{L}}(\mbox{\bf f},\mbox{\bf g},\mbox{\boldmath$\theta$})$ that is
not necessarily positive definite (Aitkin and Rubin, 1985; Quinn et al.,
1987). Nevertheless, we can use the same arguments as those used in (Keribin,
2000) 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_{\mathrm{BIC}}\to N_{*}\quad\mbox{a.s.},\quad\mbox{when}\quad
M,K\to\infty,$ (63)

provided variances
${\bf\Sigma}_{1},\ldots,{\bf\Sigma}_{N},{\bf\Gamma}_{1},\ldots,{\bf\Gamma}_{N}$
are non-degenerate and the sequence $\frac{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 10), choosing correctly the number of objects in all the
cases.

##  8 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 (Forsyth and Ponce, 2003) but practical methods are
strongly dependent on the lighting 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; Christensen et al.,
2007), 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 (Dibiase et
al., 2001). 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
$\mbox{\boldmath$s$}_{n}=(x_{n},y_{n},z_{n})^{\top}$.

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 vector $\mbox{\boldmath$f$}_{m}=(u_{m},v_{m},d_{m})^{\top}$
corresponds to a 2D image location $(u_{m},v_{m})$ and to an associated
binocular disparity $d_{m}$. Considering a projective camera model (Faugeras,
1993) it is straightforward to define an invertible function
$\mathcal{F}:\mathbb{R}^{3}\to\mathbb{R}^{3}$ that maps
$\mbox{\boldmath$s$}=(x,y,z)^{\top}$ onto
$\mbox{{\boldmath$f$}}=(u,v,d)^{\top}$:

$\mathcal{F}(\mbox{\boldmath$s$})=\left(\frac{x}{z},\frac{y}{z},\frac{1}{z}\right)^{\top}\quad\mbox{and
}\quad\mathcal{F}^{-1}(\mbox{\boldmath$f$})=\left(\frac{u}{d},\frac{v}{d},\frac{1}{d}\right)^{\top}.$
(64)

This model corresponds to a rectified camera pair (Hartley and Zisserman,
2000) and it can be easily generalized to more complex binocular geometries
(Hansard and Horaud, 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
$\mathcal{G}:\mathbb{R}^{3}\to\mathbb{R}$ maps
$\mbox{\boldmath$s$}=(x,y,z)^{\top}$ onto a 1D audio observation:

$g=\mathcal{G}(\mbox{\boldmath$s$})=\frac{1}{c}\Bigl{(}\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{M_{1}}\|-\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{M_{2}}\|\Bigr{)}.$
(65)

Here $c$ is the sound speed and $\mbox{\boldmath$s$}_{M_{1}}$ and
$\mbox{\boldmath$s$}_{M_{2}}$ 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 $\mathcal{F}$ and $\mathcal{G}$ and their derivatives are Lipschitz
continuous. We prove the following theorem:

######  Theorem 3.

The functions $\mathcal{F}$, $\mathcal{F}^{\prime}$, $\mathcal{G}$ and
$\mathcal{G}^{\prime}$ are Lipschitz continuous with constants
$L_{\mathcal{F}}=z_{\mathrm{min}}^{-1}\sqrt{3}$,
$L^{\prime}_{\mathcal{F}}=z_{\mathrm{min}}^{-2}$,
$L_{\mathcal{G}}=\|\mbox{\boldmath$s$}_{\mathrm{M_{1}}}-\mbox{\boldmath$s$}_{\mathrm{M_{2}}}\|(cR)^{-1}$
and $L^{\prime}_{\mathcal{G}}=3(cR)^{-1}$ in the domain
$\mathbb{S}=\\{|z|>z_{\mathrm{min}}>1\\}\cap\Big{\\{}\min\\{\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{1}}}\|,\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{2}}}\|\\}>R>1\Big{\\}}.$

Proof:  The derivatives of $\mathcal{F}$ and $\mathcal{G}$ are given by:

$\displaystyle\mathcal{F}^{\prime}(\mbox{\boldmath$s$})$ $\displaystyle=$
$\displaystyle\frac{1}{z}\left[\begin{array}[]{lll}1&0&-x/z\\\ 0&1&-y/z\\\
0&0&-1/z\end{array}\right]$ (69)
$\displaystyle\mathcal{G}^{\prime}(\mbox{\boldmath$s$})$ $\displaystyle=$
$\displaystyle\frac{1}{c}\left(\frac{\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{1}}}}{\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{1}}}\|}-\frac{\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{2}}}}{\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{2}}}\|}\right).$
(70)

The eigenvalues of $\mathcal{F}^{\prime}(\mbox{\boldmath$s$})$ are $1/z$ and
$-1/z^{2}$, so
$\|\mathcal{F}^{\prime}(\mbox{\boldmath$s$})\|\leq\max\\{z^{-1},z^{-2}\\}\leq
z_{\mathrm{min}}^{-1}$, from which it follows that $L_{\mathcal{F}}$ can be
taken as $L_{\mathcal{F}}=z_{\mathrm{min}}^{-1}\sqrt{3}$. Also
$\|\mathcal{F}^{\prime}(\mbox{\boldmath$s$})-\mathcal{F}^{\prime}(\tilde{\mbox{\boldmath$s$}})\|\leq\max\\{|z^{-1}-\tilde{z}^{-1}|,|z^{-2}-\tilde{z}^{-2}|\\}\leq
z_{\mathrm{min}}^{-2}\|\mbox{\boldmath$s$}-\tilde{\mbox{\boldmath$s$}}\|$, so
that $L^{\prime}_{\mathcal{F}}$ can be set to
$L^{\prime}_{\mathcal{F}}=z_{\mathrm{min}}^{-2}$.

Introducing
$\mbox{\boldmath$e$}_{1}=\frac{\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M1}}}{\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{1}}}\|}$
and
$\mbox{\boldmath$e$}_{2}=\frac{\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M2}}}{\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{2}}}\|}$,
it comes $\|\mbox{\boldmath$e$}_{1}\|=\|\mbox{\boldmath$e$}_{2}\|=1$ and
$\mathcal{G}^{\prime}(\mbox{\boldmath$s$})=\frac{1}{c}(\mbox{\boldmath$e$}_{1}-\mbox{\boldmath$e$}_{2})$.
Provided that $\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{1}}}\|$
and $\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{2}}}\|$ are both
greater than $R$, it follows
$\|\mathcal{G}^{\prime}(\mbox{\boldmath$s$})\|=\frac{1}{c}\|\mbox{\boldmath$e$}_{1}-\mbox{\boldmath$e$}_{2}\|\leq\|\mbox{\boldmath$s$}_{\mathrm{M_{1}}}-\mbox{\boldmath$s$}_{\mathrm{M_{2}}}\|(cR)^{-1}$
and so
$L_{\mathcal{G}}=\|\mbox{\boldmath$s$}_{\mathrm{M_{1}}}-\mbox{\boldmath$s$}_{\mathrm{M_{2}}}\|(cR)^{-1}$.
Then, the second derivative of $\mathcal{G}$ is given by

$\mathcal{G}^{\prime\prime}(\mbox{\boldmath$s$})=\frac{1}{c\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{1}}}\|}(\mbox{\bf
I}-\mbox{\boldmath$e$}_{1}\mbox{\boldmath$e$}_{1}^{\top})-\frac{1}{c\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{2}}}\|}(\mbox{\bf
I}-\mbox{\boldmath$e$}_{2}\mbox{\boldmath$e$}_{2}^{\top}).$

so that
$\|\mathcal{G}^{\prime\prime}(\mbox{\boldmath$s$})\|\leq\left|\frac{1}{c\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M1}}\|}-\frac{1}{c\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M2}}\|}\right|+\sup\limits_{\|\mbox{\boldmath$v$}\|=1}\frac{2\mbox{\boldmath$e$}_{1}\mbox{\boldmath$e$}_{1}^{\top}\mbox{\boldmath$v$}}{c\min\\{\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M1}}\|,\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M2}}\|\\}}\leq
3(cR)^{-1}$, and $L^{\prime}_{\mathcal{G}}$ can be set to
$L^{\prime}_{\mathcal{G}}=3(cR)^{-1}.\quad\blacksquare$

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 $\tilde{\mbox{\boldmath$s$}}^{(\nu)}$.
Equation (61) gives the largest possible step size $\rho^{(\nu)}$, so setting
$z_{\mathrm{min}}^{(\nu)}=z^{(\nu)}-\rho^{(\nu)}$ and
$R^{(\nu)}=\min\\{\|\tilde{\mbox{\boldmath$s$}}^{(\nu)}-\mbox{\boldmath$s$}_{\mathrm{M_{2}}}\|,\|\tilde{\mbox{\boldmath$s$}}^{(\nu)}-\mbox{\boldmath$s$}_{\mathrm{M_{1}}}\|\\}-\rho^{(\nu)}$,
provides local Lipschitz constants that insure the update not to quit
$\mathbb{S}^{(\nu)}=\\{|z|>z^{(\nu)}_{\mathrm{min}}\\}\cap\Big{\\{}\min\\{\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{1}}}\|,\|\mbox{\boldmath$s$}-\mbox{\boldmath$s$}_{\mathrm{M_{2}}}\|\\}>R^{(\nu)}\Big{\\}}$.
Second, we propose four possibilities to set the initial object parameter
values $\tilde{\mbox{\boldmath$s$}}^{(0)}_{n}$: (i) it can be taken to be the
previously estimated object position $\mbox{\boldmath$s$}^{(q-1)}_{n}$, (ii)
it can be set to $\mathcal{F}^{-1}(\bar{\mbox{\boldmath$f$}})$ (as soon as
$\mathcal{F}$ is injective in $\mathbb{S}$), (iii) it can be found through
sampling of the manifold $\mathcal{G}^{-1}(\bar{\mbox{\boldmath$g$}})$ by
selecting the sampled value which gives the largest $Q$ value, or (iv)
similarly through sampling directly in $\mathbb{S}$. Comparisons are reported
in the following sections.

##  9 Experiments with Simulated Data

Our algorithm is first illustrated on simulated data. For simplicity we
consider $(u,d)$ and $(x,z)$ coordinates so that
$\mathbb{F}\subseteq\mathbb{R}^{2}$ and $\mathbb{S}\subseteq\mathbb{R}^{2}$.
Notice however that this preserves the projective nature of the mapping
$\mathcal{F}$, it does not qualitatively affect the results and allows to
better understand the algorithm performance. We consider three objects defined
in $\mathbb{S}$ by $\mbox{\boldmath$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 $\mbox{{\boldmath$s$}}_{1}=(-300,1000)$,
$\mbox{{\boldmath$s$}}_{2}=(10,800)$ and
$\mbox{{\boldmath$s$}}_{3}=(500,1500)$. In the PoorSep case, the coordinates
are respectively $\mbox{{\boldmath$s$}}_{1}=(-300,1000)$,
$\mbox{{\boldmath$s$}}_{2}=(10,800)$ and
$\mbox{{\boldmath$s$}}_{3}=(100,1500)$. The data in both observation spaces
$\mathbb{F}$ and $\mathbb{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
$\mathcal{F}(\mbox{\boldmath$s$}_{n})$ and
$\mathcal{G}(\mbox{\boldmath$s$}_{n})$, $n=1,2,3$. An example of simulated
data for the three mentioned configurations is shown in Figure 2, i.e.,
$(u,d)$ locations of the visual observations and ITD values of the auditory
observations.

|  |
---|---|---
|  |
(a) GoodSep | (b) PoorSep | (c) PoorPrec
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.

#### 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
${\mathcal{L}}(\mbox{\bf f},\mbox{\bf g},\mbox{\boldmath$\theta$})$ for
initial parameters $\mbox{\boldmath$\theta$}^{(0)}_{\mathrm{OSC}}$ and
$\mbox{\boldmath$\theta$}^{(0)}_{\mathrm{PSC}}$ chosen by OSC and PSC
strategies respectively. For the total number of clusters $N=1,\ldots,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.

|  |
---|---|---
(a) GoodSep | (b) PoorSep | (c) PoorPrec
Figure 3: Means and variances of log-likelihood values
${\mathcal{L}}(\mbox{\bf f},\mbox{\bf g},\mbox{\boldmath$\theta$})$ for
initial parameters $\mbox{\boldmath$\theta$}^{(0)}_{\mathrm{OSC}}$ and
$\mbox{\boldmath$\theta$}^{(0)}_{\mathrm{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.

#### Optimization.

We compared several versions of the algorithm based on various Choose and
Local Search strategies. For the initial values
$\tilde{\mbox{\boldmath$s$}_{n}}^{(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
$\bar{\mbox{\boldmath$f$}}_{n}$, calculated in the current E-step of the
algorithm for every $n$, was mapped to the parameter space and
$\tilde{\mbox{\boldmath$s$}_{n}}^{(0)}$ set to
$\tilde{\mbox{\boldmath$s$}_{n}}^{(0)}=\mathcal{F}^{-1}(\bar{\mbox{\boldmath$f$}}_{n})$
using the injectivity of $\mathcal{F}$.

  * •

When initializing from audio data (IA),
$\mathcal{G}^{-1}(\bar{\mbox{\boldmath$g$}}_{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 (Zhigljavsky, 1991) 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 $\mathbb{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 non-
decreasing 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.

Table 1:  Acronyms used for five variants of the conjugate EM algorithm. Variants correspond to different choices for the Choose and Local search procedures. Acronym |  $\tilde{\mbox{\boldmath$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

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: well-separated object but poor observation precision.

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 $\tilde{\mbox{\boldmath$s$}}^{(0)}$ is taken from the
previous EM iteration. When $\tilde{\mbox{\boldmath$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
$\mathcal{F}^{-1}(\bar{\mbox{\boldmath$f$}}_{n})$ and
$\mathcal{G}^{-1}(\bar{\mbox{\boldmath$g$}}_{n})$ with possible small
perturbations to find the best $\tilde{\mbox{\boldmath$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).

Figure 5: IAAA algorithm: parameter evolution and assignment results for the GoodSep 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. Observations assignments are depicted by different markers ($\circ$, $*$ and $\times$ 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 6: IAAA algorithm: parameter evolution for the GoodSep case in object space. 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. |  |
---|---|---
(a) GoodSep | (b) PoorSep | (c) PoorPrec
Figure 7: BIC score graphs for the three object configurations, evaluated for
models trained for different total number of clusters $N$.

Parameter evolution trajectories for the IAAA version of the algorithm in the
GoodSep case are shown in Figures 5-6. 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.

Table 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 | $\mbox{\boldmath$s$}=(x,z)$ | $\hat{\mbox{\boldmath$s$}}=(\hat{x},\hat{z})$ | $e_{\mathrm{a}}=\|\hat{\mbox{\boldmath$s$}}-\mbox{\boldmath$s$}\|$ | $e_{\mathrm{r}}=\|\hat{\mbox{\boldmath$s$}}-\mbox{\boldmath$s$}\|/\|\mbox{\boldmath$s$}\|$
GoodSep | Object 1 |  $(-300,1000)$ |  $(-300.13,997.81)$ | 2.2 | $2.1\cdot 10^{-3}$
Object 2 |  $(10,800)$ |  $(9.28,804.46)$ | 4.52 | $5.7\cdot 10^{-3}$
Object 3 |  $(500,1500)$ |  $(513.56,1555.23)$ | 56.86 | $3.5\cdot 10^{-2}$
PoorSep | Object 1 |  $(-300,1000)$ |  $(-307.47,1028.38)$ | 29.35 | $2.8\cdot 10^{-2}$
Object 2 |  $(10,800)$ |  $(14.19,895.69)$ | 95.79 | $1.2\cdot 10^{-1}$
Object 3 |  $(100,1500)$ |  $(105.02,1447.49)$ | 52.75 | $3.5\cdot 10^{-2}$
PoorPrec | Object 1 |  $(-300,1000)$ |  $(-208.86,698.51)$ | 314.97 | $0.3$
Object 2 |  $(10,800)$ |  $(8.44,703.97)$ | 96.04 | $1.2\cdot 10^{-1}$
Object 3 |  $(500,1500)$ |  $(507.65,1533.8)$ | 34.66 | $2.2\cdot 10^{-2}$
Visual Space | $\mbox{\boldmath$f$}=(u,d)$ | $\hat{\mbox{\boldmath$f$}}=(\hat{u},\hat{d})$ | $e_{\mathrm{a}}=\|\hat{\mbox{\boldmath$f$}}-\mbox{\boldmath$f$}\|$ | $e_{\mathrm{r}}=\|\hat{\mbox{\boldmath$f$}}-\mbox{\boldmath$f$}\|/\|\mbox{\boldmath$f$}\|$
GoodSep | Object 1 |  $(-0.3,0.001)$ |  $(-0.3008,0.001)$ | $7.87\cdot 10^{-4}$ | $2.6\cdot 10^{-3}$
Object 2 |  $(0.0125,0.00125)$ |  $(0.0115,0.00124)$ | $9.59\cdot 10^{-4}$ | $7.6\cdot 10^{-2}$
Object 3 |  $(0.3333,0.00067)$ |  $(0.3302,0.00064)$ | $31.21\cdot 10^{-4}$ | $9.3\cdot 10^{-3}$
PoorSep | Object 1 |  $(-0.3,0.001)$ |  $(-0.299,0.001)$ | $1.02\cdot 10^{-3}$ | $3.4\cdot 10^{-3}$
Object 2 |  $(0.0125,0.00125)$ |  $(0.0159,0.00112)$ | $3.36\cdot 10^{-3}$ | $2.6\cdot 10^{-1}$
Object 3 |  $(0.6667,0.00067)$ |  $(0.7131,0.00238)$ | $4.95\cdot 10^{-3}$ | $7.4\cdot 10^{-2}$
PoorPrec | Object 1 |  $(-0.3,0.001)$ |  $(-0.299,0.0014)$ | $10.8\cdot 10^{-4}$ | $3.5\cdot 10^{-3}$
Object 2 |  $(0.0125,0.00125)$ |  $(0.012,0.00142)$ | $5.38\cdot 10^{-4}$ | $4.3\cdot 10^{-2}$
Object 3 |  $(0.3333,0.00067)$ |  $(0.331,0.00065)$ | $23.56\cdot 10^{-4}$ | $7.1\cdot 10^{-3}$
Audio Space | $g$ | $\hat{\mbox{\boldmath$g$}}$ | $e_{\mathrm{a}}=|\hat{\mbox{\boldmath$g$}}-\mbox{\boldmath$g$}|$ | $e_{\mathrm{r}}=|\hat{\mbox{\boldmath$g$}}-\mbox{\boldmath$g$}|/|\mbox{\boldmath$g$}|$
GoodSep | Object 1 | $-49.71$ | $-49.8$ | $0.09$ | $1.9\cdot 10^{-3}$
Object 2 | $-8.22$ | $-8.35$ | $0.13$ | $1.6\cdot 10^{-2}$
Object 3 | $34.75$ | $34.37$ | $0.38$ | $1.1\cdot 10^{-2}$
PoorSep | Object 1 | $-49.71$ | $-49.59$ | $0.12$ | $2.3\cdot 10^{-3}$
Object 2 | $-8.22$ | $-7.76$ | $0.46$ | $5.6\cdot 10^{-2}$
Object 3 | $-0.66$ | $-0.02$ | $0.65$ | $9.7\cdot 10^{-1}$
PoorPrec | Object 1 | $-49.71$ | $-49.49$ | $0.22$ | $4.4\cdot 10^{-3}$
Object 2 | $-8.22$ | $-8.28$ | $0.06$ | $7.6\cdot 10^{-3}$
Object 3 | $34.75$ | $34.47$ | $0.29$ | $8.3\cdot 10^{-3}$

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
$\hat{\mbox{\boldmath$s$}}=(\hat{x},\hat{z})$,
$\hat{\mbox{\boldmath$f$}}=(\hat{u},\hat{d})$ and $\hat{\mbox{\boldmath$g$}}$
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 10, 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, $\mathcal{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
$\mu_{\mathcal{F}}$ and variance $\Sigma_{\mathcal{F}}$, it follows that
through the linear approximation of $\mathcal{F}^{-1}$, the variance in object
space is
$\frac{\partial\mathcal{F}^{-1}(\mu_{\mathcal{F}})}{\partial\mbox{\boldmath$f$}}{\bf\Sigma}_{\mathcal{F}}\frac{\partial\mathcal{F}^{-1}(\mu_{\mathcal{F}})}{\partial\mbox{\boldmath$f$}}^{\top}$.
Then, the $z$ coordinate covariance for an object $n$ is approximately
proportional to the $d$ covariance for the object multiplied by $z_{n}^{4}$.
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)
to the models, trained for that $N$. The BIC score graphs are shown on Figure
7. The total number of objects $N$ is correctly determined in all the 3 cases
of object configurations, from which we conclude that the BIC criterion
provides reliable model selection in our case.

##  10 Experiments with Real Data

In this section we evaluate the effectiveness of our algorithms in estimating
the 3D locations of AV objects, i.e., a person localization task. The examples
used below are from a database of realistic AV scenarios described in detail
in (Arnaud et al., 2008).

The experimental setup consists of a mannequin equipped with a pair of
microphones fixed into its ears and a pair of stereoscopic cameras mounted
onto its forehead (this device was developed within the
POP222http://perception.inrialpes.fr/POP/ project). Each data set comprises
two audio tracks, two image sequences, as well as the calibration information.
All the recordings were performed in an ordinary room with no special
adjustments to its acoustics or appearance. Thus the data contain both visual
background information, and auditory noise, reverberations in particular. This
configuration best mimics what a person would hear and see in a standard
indoor environment.

We tested our multimodal clustering method with two scenarios: a meeting and a
cocktail party, Table 3:

  * •

The meeting scenario333
http://perception.inrialpes.fr/CAVA_Dataset/Site/data.html#M1 is a recording
of a discussion held by five persons sitting around a table, only three of
them being visible. It lasts 25 seconds and contains a total of about 8000
visual and 600 audio observations. The three visible persons perform head and
body movements while taking speech turns. Sometimes two persons (visible or
not) speak simultaneously.

  * •

The cocktail party scenario444
http://perception.inrialpes.fr/CAVA_Dataset/Site/data.html#CTMS3 shows a
dynamic scene with three persons walking in a room and taking speech turns.
Occasionally, one speaker is hidden by another person and two persons may
speak simultaneously. Speakers may go in and out of the two cameras field of
view. Moreover, there are sounds emitted by the persons’ steps. The recording
lasts 30 seconds and contains a total of about 12500 visual and 3400 audio
observations.

scenario |  visible persons |  speaking persons |  visual background |  audio noise |  occluded speakers |  audio overlap
---|---|---|---|---|---|---
meeting | 3 | 5 | yes | yes | no | yes
cocktail party | 3 | 3 | yes | yes | yes | yes
Table 3:  Summary of the main characteristics of the two scenarios used to
evaluate the multimodal clustering algorithm.

###  10.1 Preprocessing and Algorithm Initialization

Visual observations, $\mathbf{f}$, are obtained as follows. First we detect
points of interest (POI) in both the left and right images and we select those
points that correspond to a moving scene object. Second we perform stereo
matching such that a disparity value is associated with each matched point.

In practice we used the POI detector described in (Harris and Stephens, 1988).
This detector is known to have high repeatability in the presence of texture
and to be photometric invariant. We analyse each image point detected this way
and we select those points associated with a significant motion pattern.
Motion patterns are obtained in a straightforward manner. A temporal intensity
variance $\sigma_{t}$ is estimated at each POI. Assuming stable lighting
conditions, the POI belongs to a static scene object if its temporal intensity
variance is low and non-zero due to a camera noise only. For image points
belonging to a dynamic scene object, the local variance is higher and depends
on the texture of the moving object and on the motion speed. In our
experiments, we estimated the local temporal intensity variance $\sigma_{t}$
at each POI, from a collection of 5 consecutive frames. The point is labelled
“motion” if $\sigma_{t}>5$ (for 8-bit gray-scale images), otherwise it is
labelled as “static”. The motion-labelled points are then matched and the
associated disparities are estimated using standard stereo methods. In
practice the results shown in this paper are obtained with the method
described in (Hansard and Horaud, 2007) using the INTEL’s OpenCV camera
calibration software 555http://www.intel.com/technology/computing/opencv.
Overall, this provides the $(u,v,d)^{\top}$ to $(x,y,z)^{\top}$ mapping (64).
Examples are shown on Figure 8. Alternatively, we could have used the
spatiotemporal point detector described in (Laptev, 2005). This methods is
designed to detect points in a video stream having large local variance in
both the spatial and temporal domains, thus representing abrupt events in the
stream. However, such points are quite rare in our dataset.

(a) Meeting: There are five speakers but only three are visible. (b) Cocktail
party: The three speakers walk in the room.

Figure 8:  This figure illustrates how the audio-visual data are preprocessed.
Visual points of interest (POI) associated with scene motion are matched
between the left and right images. The histograms of the interaural time
difference (ITD) observations correspond to a “segment” of 0.3 seconds. The
audio-visual calibration allows us to filter out auditory data that falls
outside the field of view of the two cameras. Notice the large number of
auditory perturbations corresponding to noise, reverberations, as well as to
speakers that are outside the visual field of view. In these examples, there
are two simultaneous speakers: (a) S1 and S4 and (b) S1 and S2. Notice that S4
is easily eliminated because its associated ITD falls outside the visual field
of view.

Auditory observations, $\mathbf{g}$, are obtained as follows. Our method uses
interaural time differences (ITD) which are detected through the analysis of
the cross-correlogram of the filtered left- and right-microphone one-
dimensional signals for every frequency band (Christensen et al., 2007). Like
any other audio-visual fusion method, one needs to perform audio-visual
calibration, namely to estimate the positions of the microphones and the
positions and orientations of the cameras in a common world coordinate system.
This is done using the method described in (Arnaud et al., 2008).

In order to initialize the algorithm’s parameter values we used the Parameter
Space Candidates (PSC) initialization strategy described in Section 6.
Although real-data distributions do not strictly correspond to the case of
Gaussian mixtures, the initialization strategy that we have adopted remains
relevant. This originates from the fact that parameter space sampling with
configuration restrictions plays the role of a global optimization method
similar to Monte-Carlo sampling in the method of generations (Zhigljavsky and
Žilinskas, 2008). It helps to avoid local maxima and allows to quickly find a
set of appropriate initial parameters. Local distribution density modes occur
to be good candidates to initialize cluster centers. As in the case of
simulated data, we used the BIC score, i.e., Section 7 to select the optimal
number of audio-visual clusters.

###  10.2 Results and Discussion

The experimental validation described below was performed with two goals in
mind. Firstly, we wanted to check that our method was stable and robust with
real data gathered in complex situations, that it correctly finds the number
of clusters and that it efficiently determines the model’s parameters, i.e.,
the 3D positions of the audio-visual objects composing a scene. Secondly, we
wanted to test the model’s capability to deal with dynamic changes in the
scene, yet in the presence of acoustic noise/reverberations and visually
occluded persons, etc. Below we provide a detailed account of the results
obtained with the meeting and cocktail-party audio-visual sequences.

The audio-visual recordings are split into “segments”, each segment lasts 0.3
seconds. At 25 frames/second this corresponds to approximately eight video
frames. The initialization method described in Section 6 and the model
selection method described in Section 7 are combined and applied to the first
segment in order to find initial parameter values and to estimate the number
of components (the number of audio-visual objects) to be used by the conjugate
EM algorithm. Consequently, the parameters estimated for one segment are used
to initialize the parameters for the next segment, while the number of
components remains constant.

  * •

Quasi-static scene. The meeting situation corresponds to the well-separated
case which is referred to as GoodSep in the previous section. The
initialization strategy performs well and the candidate configuration obtained
by the initialization step is relatively close to the optimal one found by the
EM algorithm described in detail in Section 3.3. In fact, the likelihood
evolution reported in Figure 9 shows that convergence is reached in about 20
iterations of EM, which is three times faster than in the simulated GoodSep
case reported in Figure 4. The 3D position estimates are quite accurate, in
particular the natural alignment of the speakers along the table is clearly
seen in the $XZ$ plane. Even though in practice, the data are not piecewise
Gaussian and the outliers are not uniformly distributed, our method performs
quite well, which illustrates its robustness when dealing with real-data
distributions. Figure 10 shows sequential results obtained in this case. The
speech sources are correctly detected even in the case when two persons are
simultaneously active.

  * •

Dynamic scene. The cocktail party situation corresponds to the partially
occluded case which is referred to as PoorSep in the previous section. In this
case, the locations of the audio-visual objects varies over time, as well as
their number. Nevertheless, we assume that these changes are rather slow. We
did not attempt to tune our algorithm to the dynamic case. Hence, we use the
same initialization strategy as in the quasi-static case which is briefly
summarized above. Figure 11 shows the results obtained in this case.

Figure 9: An example of applying the proposed EM algorithm to a time interval of 20 seconds of the meeting scenario. The results are shown in the visual and auditory observation spaces as well as in the parameter space. The initial parameter values are shown with three stars while the parameter evolution trajectories are shown with circles of increasing size. The final observation-to-cluster assignments are shown in color: red, blue, and green for the three Gaussian components and light-blue for the outlier component. The log-likelihood curve (bottom-right) shows that the algorithm converged after 20 iterations.  |
---|---
(a) frames 1001-1010 | (b) frames 1011-1020
|
(c) frames 1021-1030 | (d) frames 1031-1040
|
(e) frames 1041-1050 | (f) frames 1051-1060
Figure 10:  Results obtained in the case of the meeting scenario shown overlapped onto the left image. Sixty frames (1001 to 1060) were split into six segments. Parameter initialization and model selection were performed on the first segment (frames 1-10) and are not shown. The “visual” covariance matrices associated with the 3 Gaussian components are projected onto the image plane. The white dots correspond to the projected 3D locations estimated by the algorithm. The blue, green, and red colors encode the observation-to-cluster assignments and the active speaker is marked with a corresponding symbol. The algorithm correctly estimates speech sources, even in the case when two speakers are active.  |
---|---
(a) frames 181-190 | (b) frames 191-200
|
(c) frames 201-210 | (d) frames 211-220
|
(e) frames 221-230 | (f) frames 231-240
Figure 11:  Results obtained in the case of the cocktail party scenario shown
overlapped onto the left image. As in the previous case, sixty frames (181 to
240) were split into six segments. Parameter initialization and model
selection were performed on the first segment (frames 1-10) and are not shown.
As expected, well separated objects, (a)-(c), are correctly handled. While
partial occlusion, (d)-(e) is also handled correctly, the algorithm fails to
deal with a complete occlusion, (f).

Overall, the proposed method performs well on data collected in a natural
environment. The initialization strategy and the model selection criterion
proved to be robust to noise and to minor deviations from the Gaussian
distribution assumption. It possesses the features of a global optimization
method which enables to find initial parameter values that are close to
optimal ones. In both examples, the parameter initialization and model
selection were performed on the first audio-visual data segment. This
certainly biases the overall results. Indeed, in both cases, the
initialization and model selection algorithms dealt with a case were the
objects were well separated. One could rerun initialization and model
selection on every data segment, at the cost of a less efficient procedure.

The conjugate clustering method automatically weights the auditory and visual
modalities, in terms of precision and amount of observations, to infer the
parameter values. We noticed that, in general, the visual data are considered
by the algorithm as more reliable. This can be explained by the fact that, in
practice, the auditory signals are contaminated with noise and reverberations.
This typically smooths the histogram peaks in the ITD domain and adds false
peaks, as can be seen in Figures 10 and 11. As reverberations are natural for
most of the environments and sound sources, we added auditory cluster
variances to model the local smoothing effect, as well as an outlier category
to treat false peaks. In general, if the data is gathered using a small time
interval, reverberations and noise have higher effect, the observations are
scattered and auditory spatial localization is poor. At the same time,
widening the time interval would result in sharper peaks for sound sources
that are smoothed due to reverberations and dynamics of the scene, and hence
the auditory temporal localization will be less accurate. Thus the auditory
data are typically sparse both in time and space. The temporal discontinuity
of the auditory data together with the lack of resolution makes it less
reliable than the visual data.

Although our multimodal clustering model has no built-in dynamic capability,
as is the case with target-tracking methods based on the Kalman filter, the
implemented algorithm performs quite well in the case of partial visual
occlusions, as illustrated in the cocktail party scenario.

##  11 Conclusions

We proposed a novel framework 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 $\mathbb{F}_{1},\ldots,\mathbb{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 $\mathbb{F}$ and
$\mathbb{G}$, would have to be replaced by a sum of $J$ terms corresponding to
$\mathbb{F}_{1},\ldots,\mathbb{F}_{J}$ in the formulas of Section 3.

In particular, adding Gaussian priors on parameters (i.e., priors, covariance
matrices and objet locations) would not essentially change the formulae. For a
large class of dynamics equations, the update expressions (29)-(32) for priors
and variances will remain in closed form, whereas the function
$Q_{n}^{{(q)}}(\mbox{\boldmath$s$})$ in (34) will receive an additional term
$\log P(\mbox{\boldmath$s$})$. For instance, multimodal dynamic inference of
parameter values for Brownian dynamics (van Kampen, 2007) can be performed by
means of the formulated model. Gaussian priors would add a quadratic term
similar to the others in (34), that can be viewed as an ‘observation’ from the
ambient space modality. Thus the optimization algorithm would not require any
changes and would give an unbiased estimate.

Also, the assumption that assignment variables $a$ and $b$ are independent
could be relaxed. An appropriate approach to perform inference in a non
independent case would be to consider variational approximations (Jordan et
al., 1998) and in particular a variational EM (VEM) framework. The general
idea would be to approximate the joint distribution $P(\mbox{\boldmath$a$})$
by a distribution from a restricted class of probability distributions that
factorize as
$\tilde{P}(\mbox{\boldmath$a$})=\prod\limits_{m=1}^{M}\tilde{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 modelling 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
heterogenous information sources. Therefore, it has advantages over other
methods that include ad-hoc processing while being open to incorporation of
more task dependent information.

## 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.

## References

  * Aitkin and Rubin (1985) Aitkin, M. and Rubin, D. (1985).  Estimation and hypothesis testing in finite mixture models.  Journal of the Royal Statistical Society. Series B (Methodological), 47(1), 67–75.
  * Allen (1995) Allen, P. (1995).  Integrating vision and touch for object recognition tasks.  In Luo and Kay, editors, Multisensor Integration and Fusion for Intelligent Machines and Systems, pages 407–440. Ablex Publishing Corporation.
  * Anastasio et al. (2000) Anastasio, T. J., Patton, P. E., and Belkacem-Boussaid, K. E. (2000).  Using Bayes’ rule to model multisensory enhancement in the superior colliculus.  Neural Computation, 12(5), 1165–1187.
  * Arnaud et al. (2008) Arnaud, E., Christensen, H., Lu, Y., Barker, J., Khalidov, V., Hansard, M., Holveck, B., Mathieu, H., Narasimha, R., Forbes, F., and Horaud, R. (2008).  The CAVA corpus: Synchronized stereoscopic and binaural datasets with head movements.  In Proc. of ACM/IEEE Tenth International Conference on Multimodal Interfaces.
  * Beal et al. (2003) Beal, M., Jojic, N., and Attias, H. (2003).  A graphical model for audiovisual object tracking.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 25(7), 828–836.
  * Bishop (2006) Bishop, C. (2006).  Pattern Recognition and Machine Learning.  Springer.
  * Boyles (1983) Boyles, R. (1983).  On the convergence of EM algorithms.  Journal of the Royal Statistical Society: Series B, 45(1), 47–50.
  * Castellanos and Tardos (1999) Castellanos, J. and Tardos, J. (1999).  Mobile Robot Localization and Map Building: A Multisensor Fusion Approach.  Kluwer, Boston, MA.
  * Celeux et al. (2003) Celeux, G., Forbes, F., and Peyrard, N. (2003).  EM procedures using mean field-like approximations for Markov model-based image segmentation.  Pattern Recognition, 36(1), 131–144.
  * Checka et al. (2004) Checka, N., Wilson, K., Siracusa, M., and Darrell, T. (2004).  Multiple person and speaker activity tracking with a particle filter.  In Proc. of IEEE Conference on Acoustics, Speech, and Signal Processing, pages 881–884. IEEE.
  * Chen and Rui (2004) Chen, Y. and Rui, Y. (2004).  Real-time speaker tracking using particle filter sensor fusion.  Proceedings of IEEE, 92(3), 485–494.
  * Christensen et al. (2007) Christensen, H., Ma, N., Wrigley, S., and Barker, J. (2007).  Integrating pitch and localisation cues at a speech fragment level.  In Proc. of Interspeech, pages 2769–2772.
  * Coiras et al. (2007) Coiras, E., Baralli, F., and Evans, B. (2007).  Rigid data association for shallow water surveys.  IET Radar, Sonar and Navigation, 1(5), 354–361.
  * Comaniciu and Meer (2002) Comaniciu, D. and Meer, P. (2002).  Mean shift: a robust approach towards feature space analysis.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 24(5), 603–619.
  * Dempster et al. (1977) Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977).  Maximum likelihood from incomplete data via the EM algorithm (with discussion).  Journal of the Royal Statistical Society: Series B, 39(1), 1–38.
  * Dibiase et al. (2001) Dibiase, J., Silverman, H., and Brandstein, M. (2001).  Robust localization in reverberant rooms.  In M. Brandstein and D. Ward, editors, Microphone Arrays: Signal Processing Techniques and Applications, chapter 8. Springer.
  * Ernst and Banks (2002) Ernst, M. O. and Banks, M. S. (2002).  Humans integrate visual and haptic information in a statistically optimal fashion.  Nature, 415, 429–433.
  * Faugeras (1993) Faugeras, O. D. (1993).  Three Dimensional Computer Vision: A Geometric Viewpoint.  MIT Press, Boston.
  * Fisher III and Darrell (2004) Fisher III, J. W. and Darrell, T. (2004).  Speaker association with signal-level audiovisual fusion.  IEEE Transactions on Multimedia, 6(3), 406–413.
  * Fisher III et al. (2001) Fisher III, J. W., Darrell, T., Freeman, W. T., and Viola, P. (2001).  Learning joint statistical models for audio-visual fusion segregation.  In Proceedings of Annual Conference on Advances in Neural Information Processing Systems, Vancouver, BC, Canada.
  * Forsyth and Ponce (2003) Forsyth, D. and Ponce, J. (2003).  Computer Vision – A Modern Approach.  Prentice Hall, New Jersey.
  * Gatica-Perez et al. (2007) Gatica-Perez, D., Lathoud, G., Odobez, J.-M., and McCowan, I. (2007).  Audiovisual probabilistic tracking of multiple speakers in meetings.  IEEE Transactions on Audio, Speech, and Language Processing, 15(2), 601–616.
  * Hall and McMullen (2004) Hall, D. L. and McMullen, S. A. H. (2004).  Mathematical Techniques in Multisensor Data Fusion.  Artech House, New York.
  * Hansard and Horaud (2007) Hansard, M. and Horaud, R. (2007).  Patterns of binocular disparity for a fixating observer.  In Proc. of Second International Symposium of Advances in Brain, Vision, and Artificial Intelligence, pages 308–317. Springer.
  * Hansard and Horaud (2008) Hansard, M. and Horaud, R. (2008).  Cyclopean geometry of binocular vision.  Journal of the Optical Society of America A, 25(9), 2357 2369.
  * Harris and Stephens (1988) Harris, C. and Stephens, M. (1988).  A combined corner and edge detector.  In Proc. of Fourth Alvey Vision Conference, pages 147–151.
  * Hartley and Zisserman (2000) Hartley, R. and Zisserman, A. (2000).  Multiple View Geometry in Computer Vision.  Cambridge University Press, Cambridge, UK.
  * Haykin and Chen (2005) Haykin, S. and Chen, Z. (2005).  The cocktail party problem.  Neural Computation, 17, 1875–1902.
  * Heckmann et al. (2002) Heckmann, M., Berthommier, F., and Kroschel, K. (2002).  Noise adaptive stream weighting in audio-visual speech recognition.  EURASIP Journal on Applied Signal Processing, 11, 1260–1273.
  * Hospedales and Vijayakumar (2008) Hospedales, T. and Vijayakumar, S. (2008).  Structure inference for Bayesian multisensory scene understanding.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 30(12), 2140–2157.
  * Jordan et al. (1998) Jordan, M., Ghahramani, Z., Jaakkola, T., and Saul, L. (1998).  An introduction to variational methods for graphical models.  In M. Jordan, editor, Learning in Graphical Models, pages 105–162. Kluwer Academic.
  * Joshi and Sanderson (1999) Joshi, R. and Sanderson, A. (1999).  Multisensor Fusion: A Minimal Representation Framework.  World Scientific.
  * Keribin (2000) Keribin, C. (2000).  Consistent estimation of the order of mixture models.  Sankhya: the Indian Journal of Statistics, Series A, 62(1), 49–66.
  * King (2004) King, A. J. (2004).  The superior colliculus.  Current Biology, 14, 335–338.
  * King (2005) King, A. J. (2005).  Multisensory integration: strategies for synchronization.  Current Biology, 15, 339–341.
  * Kushal et al. (2006) Kushal, A., Rahurkar, M., Fei-Fei, L., Ponce, J., and Huang, T. (2006).  Audio-visual speaker localization using graphical models.  In Proc. of the Eighteenth International Conference on Pattern Recognition, pages 291–294.
  * Laptev (2005) Laptev, I. (2005).  On space-time interest points.  Int. J. of Comp. Vis., 64(2-3), 107–123.
  * Majumder et al. (2001) Majumder, S., Scheding, S., and Durrant-Whyte, H. (2001).  Multisensor data fusion for underwater navigation.  Robotics and Autonomous Systems, 35, 97–108.
  * McLachlan and Krishnan (1996) McLachlan, G. and Krishnan, T. (1996).  The EM algorithm and extensions.  Wiley.
  * McLachlan and Peel (2000) McLachlan, G. J. and Peel, D. (2000).  Finite Mixture Models.  Wiley, New-York.
  * Mitchell (2007) Mitchell, H. (2007).  Multi-Sensor Data Fusion.  Springer, Berlin Heidelberg.
  * Naus and van Wijk (2004) Naus, H. and van Wijk, C. (2004).  Simultaneous localization of multiple emitters.  IEE Proceedings Radar Sonar and Navigation, 151, 65–70.
  * Nefian et al. (2002) Nefian, A. V., Liang, L., Pi, X., Liu, X., and Murphy, K. (2002).  Dynamic Bayesian networks for audio-visual speech recognition.  EURASIP Journal of Applied Signal Processing, 2002(1), 1274–1288.
  * Perez et al. (2004) Perez, P., Vermaak, J., and Blake, A. (2004).  Data fusion for visual tracking with particles.  Proceedings of IEEE, 92(3), 495–513.
  * Polyak (1987) Polyak, B. (1987).  Introduction to Optimization.  New York: Optimization Software, Publications Division.
  * Pouget et al. (2002) Pouget, A., Deneve, S., and J.-R., D. (2002).  A computational perspective on the neural basis of multisensory spatial representations.  Nature Review Neuroscience, 3, 741–747.
  * Quinn et al. (1987) Quinn, B., McLachlan, G., and N.L.Hjort (1987).  A note on the Aitkin-Rubin approach to hypothesis testing in mixture models.  Journal of the Royal Statistical Society. Series B (Methodological), 49(3), 311–314.
  * Schwarz (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  Annals of Statistics, 6, 461–464.
  * Shao and Barker (2008) Shao, X. and Barker, J. (2008).  Stream weight estimation for multistream audio-visual speech recognition in a multispeaker environment.  Speech Communication, 50(4), 337–353.
  * Smith and Singh (2006) Smith, D. and Singh, S. (2006).  Approaches to multisensor data fusion in target tracking: A survey.  IEEE Transactions on Knowledge and Data Engineering, 18, 1041–4347.
  * van Kampen (2007) van Kampen, N. (2007).  Stochastic Processes in Physics and Chemistry.  North Holland, 3 edition.
  * Wang and Brown (2006) Wang, D. and Brown, G. J., editors (2006).  Computational Auditory Scene Analysis: Principles, Algorithms, and Applications.  Wiley-IEEE Press.
  * Zhigljavsky (1991) Zhigljavsky, A. (1991).  Theory of Global Random Search.  Kluwer Academic Publishers.
  * Zhigljavsky and Žilinskas (2008) Zhigljavsky, A. and Žilinskas, A. (2008).  Stochastic Global Optimization.  Springer.