Local Component Analysis
Authors: Nicolas Le Roux, Francis Bach
ICLR 2013

Abstract: Kernel density estimation, a.k.a. Parzen windows, is a popular density estimation method, which can be used for outlier detection or clustering. With multivariate data, its performance is heavily reliant on the metric used within the kernel. Most earlier work has focused on learning only the bandwidth of the kernel (i.e., a scalar multiplicative factor). In this paper, we propose to learn a full Euclidean metric through an expectation-minimization (EM) procedure, which can be seen as an unsupervised counterpart to neighbourhood component analysis (NCA). In order to avoid overfitting with a fully nonparametric density estimator in high dimensions, we also consider a semi-parametric Gaussian-Parzen density model, where some of the variables are modelled through a jointly Gaussian density, while others are modelled through Parzen windows. For these two models, EM leads to simple closed-form updates based on matrix inversions and eigenvalue decompositions. We show empirically that our method leads to density estimators with higher test-likelihoods than natural competing methods, and that the metrics may be used within most unsupervised learning techniques that rely on such metrics, such as spectral clustering or manifold learning methods. Finally, we present a stochastic approximation scheme which allows for the use of this method in a large-scale setting.

arXiv:1109.0093v4  [cs.LG]  10 Dec 2012
Local Component Analysis
Nicolas Le Roux
nicolas@le-roux.name
Francis Bach
francis.bach@ens.fr
INRIA - SIERRA Project - Team
Laboratoire d’Informatique de l’ ´Ecole Normale Sup´ erieure
Paris, France
Abstract
Kernel density estimation, a.k.a. Parzen windows, is
a popular density estimation method, which can be
used for outlier detection or clustering. With mul-
tivariate data, its performance is heavily reliant on
the metric used within the kernel. Most earlier work
has focused on learning only the bandwidth of the
kernel (i.e., a scalar multiplicative factor). In this
paper, we propose to learn a full Euclidean metric
through an expectation-minimisation (EM) procedure,
which can be seen as an unsupervised counterpart to
neighbourhood component analysis (NCA). In order
to avoid overﬁtting with a fully nonparametric den-
sity estimator in high dimensions, we also consider
a semi-parametric Gaussian-Parzen density model,
where some of the variables are modelled through a
jointly Gaussian density, while others are modelled
through Parzen windows. For these two models, EM
leads to simple closed-form updates based on matrix
inversions and eigenvalue decompositions. We show
empirically that our method leads to density estima-
tors with higher test-likelihoods than natural compet-
ing methods, and that the metrics may be used within
most unsupervised learning techniques that rely on lo-
cal distances, such as spectral clustering or manifold
learning methods. Finally, we present a stochastic ap-
proximation scheme which allows for the use of this
method in a large-scale setting.
1 Introduction
Most unsupervised learning methods rely on a metric
on the space of observations. The quality of the met-
ric directly impacts the performance of such techniques
and a signiﬁcant amount of work has been dedicated to
learning this metric from data when some supervised
information is available [27, 16, 2]. However, in a fully
unsupervised scenario, most practitioners use the Ma-
halanobis distance obtained from principal component
analysis (PCA). This is an unsatisfactory solution as
PCA is essentially a global linear dimension reduction
method, while most unsupervised learning techniques,
such as spectral clustering or manifold learning, are
local.
In this paper, we cast the unsupervised metric learning
as a density estimation problem with a Parzen win-
dows estimator based on a Euclidean metric. Using
the maximum likelihood framework, we derive in Sec-
tion 3 an expectation-minimisation (EM) procedure
that maximizes the leave-one-out log-likelihood, which
may be considered as an unsupervised counterpart to
neighbourhood component analysis (NCA) [16]. As
opposed to PCA, which performs a whitening of the
data based on global information, our new algorithm
globally performs a whitening of the data using only
local information, hence the denomination local com-
ponent analysis (LCA).
Like all non-parametric density estimators, Parzen
windows density estimation is known to overﬁt in high
dimensions [25], and thus LCA should also overﬁt. In
order to keep the modelling ﬂexibility of our density es-
timator while avoiding overﬁtting, we propose a semi-
parametric Parzen-Gaussian model; following [4], we
linearly transform then split our variables in two parts,
one which is modelled through a Parzen windows es-
timator (where we assume the interesting part of the
data lies), and one which is modelled as a multivari-
ate Gaussian (where we assume the noise lies). Again,
in Section 4, an EM procedure for estimating the lin-
ear transform may be naturally derived and leads to
simple closed-form updates based on matrix inversions
and eigenvalue decompositions. This procedure con-
tains no hyperparameters, all the parameters being
learnt from data.
Since the EM formulation of LCA scales quadratically
in the number of datapoints, making it impractical
for large datasets, we introduce in Section 5 both a
stochastic approximation and a subsampling technique
1
allowing us to achieve a linear cost and thus to scale
LCA to much larger datasets.
Finally, in Section 6, we show empirically that our
method leads to density estimators with higher test-
likelihoods than natural competing methods, and that
the metrics may be used within unsupervised learn-
ing techniques that rely on such metrics, like spectral
clustering.
2 Previous work
Many authors aimed at learning a Mahalanobis dis-
tance suited for local learning. While some techniques
required the presence of labelled data [16, 27, 2], oth-
ers proposed ways to learn the metric in a purely un-
supervised way, e.g., [28] who used the distance to the
k-th nearest neighbour as the local scaling around each
datapoint. Most of the other attempts at unsuper-
vised metric learning were developed in the context
of kernel density estimation, a.k.a. Parzen windows.
The Parzen windows estimator [21] is a nonparametric
density estimation model which, given n datapoints
{x1, . . . , xn}in Rd, deﬁnes a mixture model of the
form p(x) = 1
n
∑ n
j=1 K(x, xj , θ ) where K is a kernel
with compact support and parameters θ. We relax the
compact support assumption and choose K to be the
normal kernel, that is
p(x) = 1
n
n∑
j=1
N(x, xj , Σ)
∝ 1
n
√
|Σ |
n∑
j=1
exp
[
−1
2(x −xj )⊤Σ −1(x −xj )
]
,
where Σ is the covariance matrix of each Gaussian.
As the performance of the Parzen windows estimator
is more reliant on the covariance matrix than on the
kernel, there has been a large body of work, originating
from the statistics literature, attempting to learn this
matrix. However, almost all attempts are focused on
the asymptotic optimality of the estimators obtained
with little consideration for the practicality in high
dimensions. Thus, the vast majority of the work is
limited to isotropic matrices, reducing the problem to
ﬁnding a single scalar h [22, 13, 23, 9, 7, 20, 24], the
bandwidth, and the few extensions to the non-isotropic
cases are numerically expensive [14, 18].
An exception is the approach proposed in [26], which
is very similar to our method, as the authors learn the
covariance matrix of the Parzen windows estimator us-
ing local neighbourhoods. However, their algorithm
does not minimize a well-deﬁned cost function, mak-
ing it unsuitable for kernels other than the Gaussian
one, and the locality used to compute the covariance
matrix depends on parameters which must be hand-
tuned or cross-validated. Also, the modelling of all
the dimensions using the Parzen windows estimator
makes the algorithm unsuitable when the data lie on
a high-dimensional manifold. In an extension to [26],
[3] uses a neural network to compute the leading eigen-
vectors of the local covariance matrix at every point
in the space, then uses these matrices to do density
estimation and classiﬁcation. Despite the algorithm’s
impressive performance, it does not correspond to a
linear reparametrisation of the space and thus cannot
be used as a preprocessing step.
3 Local Component Analysis
Seeing the density as a mixture of Gaussians, one can
easily optimize the covariances using the EM algo-
rithm [12]. However, maximizing the standard log-
likelihood of the data would trivially lead to the degen-
erate solution where Σ goes to 0 to yield a sum of Dirac
distributions. One solution to that problem is to pe-
nalize some norm of the precision matrix to prevent it
from going to inﬁnity. Another, more compelling, way
is to optimize the leave-one-out log-likelihood, where
the probability of each datapoint xi is computed un-
der the distribution obtained when xi has been re-
moved from the training set. This technique is not
new and has already been explored both in the su-
pervised [16, 15] and in the unsupervised setting [13].
However, in the latter case, the cross-validation was
then done by hand, which explains why only one band-
width parameter could be optimized 1. We will thus
use the following criterion:
L(Σ) = −
n∑
i=1
log
[ 1
n −1
∑
j̸=i
N(xi, x j , Σ)
]
(1)
≤cst −
n∑
i=1
∑
j̸=i
λij log N(xi, x j , Σ)
+
n∑
i=1
∑
j̸=i
λij log λij , (2)
with the constraints ∀i , ∑
j̸=i λij = 1. This varia-
tional bound is obtained using Jensen’s inequality.
The EM algorithm optimizes the right-hand side of
Eq. (2) by alternating between the optimisations of λ
and Σ in turn. The algorithm is guaranteed to con-
verge, and does so to a stationary point of the true
1Most of the literature on estimating the covariance matrix
discards the log-likelihood cost because of its sensitivit y to out-
liers and prefers AMISE (see, e.g., [14]). However, in all ou r
experiments, the number of datapoints was large enough so th at
LCA did not suﬀer from the presence of outliers.
2
function over Σ deﬁned in Eq. (1). At each step, the
optimal solutions are:
λ∗
ij = N(xi, x j , Σ)
∑
k̸=i N(xi, x k, Σ) if j ̸= i (3)
λ∗
ii = 0 (4)
Σ ∗ =
∑
ij λij (xi −xj )(xi −xj )T
n . (5)
The “responsibilities” λ∗
ij deﬁne the relative proximity
of xj to xi (compared to the proximity of all the xk’s
to xi) and Σ ∗ is the average of all the local covariance
matrices.
This algorithm, which we coin LCA, for local compo-
nent analysis, transforms the data to make it locally
isotropic, as opposed to PCA which makes it glob-
ally isotropic. Fig. (1) shows a comparison of PCA
and LCA on the word sequence “To be or not to be”.
Whereas PCA is highly sensitive to the length of the
text, LCA is only aﬀected by the local shapes, thus
providing a much less distorted result. 2
First, one may note at this point that Manifold Parzen
Windows [26] is equivalent to LCA with only one step
of EM. This makes Manifold Parzen Windows more
sensitive to the choice of the original covariance ma-
trix whose parameters must be carefully chosen. As we
shall see later in the experiments, running EM to con-
vergence is important to get good accuracy when using
spectral clustering on the transformed data. Second, it
is also worth noting that, similarly to Manifold Parzen
Windows, LCA can straightforwardly be extended to
the cases where each datapoint uses its own local co-
variance matrix (possibly with a smoothing term), or
where the covariance Σ ∗ is the sum of a low-rank ma-
trix and some scalar multiplied by the identity matrix.
Not only may LCA be used to learn a linear trans-
formation of the space, but it also deﬁnes a density
model. However, there are two potential negative as-
pects associated with this method. First, in high di-
mensions, Parzen windows is prone to overﬁtting and
must be regularized [25]. Second, if there are some
directions containing a small Gaussian noise, the local
isotropy will blow them up, swamping the data with
clutter. This is common to all the techniques which
renormalise the data by the inverse of some variance.
A solution to both of these issues is to consider a prod-
uct of two densities: one is a low-dimensional Parzen
windows estimator, which will model the interesting
signal, and the other is a Gaussian, which will model
the noise.
2Since both methods are insensitive to any linear
reparametrisation of the data, we do not include the origina l
data in the ﬁgure.
4 LCA with a multiplicative
Gaussian component
We now assume that there are irrelevant dimensions
in our data which can be modelled by a Gaussian. In
other words, we consider an invertible linear transfor-
mation ( BG, B L)⊤x of the data, modelling B⊤
G x as a
multivariate Gaussian and B⊤
L x through kernel density
estimation, the two parts being independent, leading
to p(x) ∝p(B⊤
G x, B ⊤
L x) = pG(B⊤
G x)pL(B⊤
L x), where
pG is a Gaussian and pL is the Parzen windows esti-
mator, i.e.,
p(xi) ∝
⏐
⏐BGB⊤
G + BLB⊤
L
⏐
⏐
1
2
n −1
×exp
[
−1
2(xi −µ)⊤BGB⊤
G (xi −µ)
]
×

∑
j̸=i
exp
[
−1
2(xi −xj )⊤BLB⊤
L (xi −xj )
] 
 ,
with ( BG, B L) a full-rank square matrix. Using EM,
we can upper-bound the negative log-likelihood:
−2
∑
i
log p(xi) ≤tr(B⊤
G CGBG) + tr( B⊤
L CLBL)
−log |BGB⊤
G + BLB⊤
L |, (6)
with
CG = 1
n
∑
i
(xi −µ)(xi −µ)⊤ ,
CL = 1
n
∑
ij
λij (xi −xj )(xi −xj )⊤ .
The matrices BG and BL minimizing the right-hand
side of Eq. (6) may be found using the following propo-
sition (see proof in the appendix):
Proposition 1 Let BG ∈ Rd×d1 and BL ∈ Rd×d2 ,
with d = d1 + d2 and B = ( BG, B L) ∈Rd×d invertible.
Consider two symmetric positive matrices M1 and M2
in Rd×d. The problem
min
BG,B L
tr B⊤
G M1BG+tr B⊤
L M2BL−log det(BGB⊤
G +BLB⊤
L )
(7)
has a ﬁnite solution only if M1 and M2 are invertible
and, if these conditions are met, reaches its optimum
at
BG = M−1/ 2
1 U+ , B L = M−1/ 2
1 U−D−1/ 2
− ,
where U+ are the eigenvectors of M−1/ 2
1 M2M−1/ 2
1 as-
sociated with eigenvalues greater than or equal to 1,
U− are the eigenvectors of M−1/ 2
1 M2M−1/ 2
1 associated
3
Figure 1: Results obtained when transforming “To be or not to be” u sing PCA (left) and LCA (right). Left:
To make the data globally isotropic, PCA awkwardly compresses the le tters horizontally. Right: Since LCA is
insensitive to the spacing between letters and to the length of the t ext, there is no horizontal compression.
with eigenvalues smaller than 1 and D− is the diagonal
matrix containing the eigenvalues of M−1/ 2
1 M2M−1/ 2
1
smaller than 1.
The resulting procedure is described in Algorithm 1,
where all dimensions are initially modelled by the
Parzen windows estimator, which empirically yielded
the best results.
Algorithm 1 LCA - Gauss
Input: X (dataset), iterMax (maximum number of
iterations), ν (regularisation)
Output: BG (Gaussian part transformation), BL
(Parzen windows transformation)
CG ← cov(X) + νId {Initialize C to the global
covariance}
IG ←C
− 1
2
G
BG = 0 , B L = chol( C−1
G ) {Assign all dimensions to
the Parzen windows estimator }
for iter = 1:iterMax do
Mij ←exp
[
−(xi−xj )⊤B⊤
L BL(xi−xj )
2
]
, Mii ←0
λij ← Mij∑
k Mik
CL ←
∑
ij λ ij (xi−xj )(xi−xj )⊤
n + νId
[V, D ] ← eig(IGCLIG) {Eigendecomposition of
IGCLIG}
t1 = max z D(z, z ) ≤1 {Cut-oﬀ between eigenval-
ues smaller and larger than 1 }
t+ = {t|t1 ≤t ≤d}, t− = {t|1 ≤t < t 1}
V+ ← V (:, t +), V− ← V (:, t −), D− =
D(t−, t −)
BL = IGV−D−1/ 2
− , BG = IGV+
end for
Relationship with ICA. Independent component
analysis (ICA) can be seen as a density model where
x = As and s has independent components (see,
e.g., [17]). In the Parzen windows framework, this cor-
responds to modelling the density of s by a product of
univariate kernel density estimators [6]. This however
causes two problems: ﬁrst, while this assumption is
appropriate in settings such as source separation, it is
violated in most settings, and having a multivariate
kernel density estimation is preferable. Second, most
algorithms are dedicated to ﬁnding independent com-
ponents which are non-Gaussian. In the presence of
more than one Gaussian dimension, most ICA frame-
works become unidentiﬁable, while our explicit mod-
elling of such Gaussian components allows us to tackle
this situation (a detailed analysis of the identiﬁability
of our Parzen/Gaussian model is out of the scope of
this paper).
Relationship with NGCA. NGCA [4] makes an
assumption similar to ours (they rather assume an
additive Gaussian noise on top of a low-dimensional
non-Gaussian signal) but uses a projection pursuit
algorithm to iteratively ﬁnd the directions of non-
Gaussianity. Unlike in FastICA, the contrast functions
used to ﬁnd the interesting directions can be diﬀerent
for each direction. However, like all projection pur-
suit algorithms, the identiﬁcation of interesting direc-
tions gets much harder in higher dimensions, as most
of them will be almost Gaussian. Our use of a non-
parametric density estimator with a log-likelihood cost
allows us to globally optimize all directions simultane-
ously and does not rely on the model being correct.
Finally, LCA estimates all its parameters from data
as opposed to NGCA which requires the number of
non-Gaussian directions to be set.
Escaping local optima. Though our model allows
for the modiﬁcation of the number of dimensions mod-
elled by the Gaussian through the analysis of the spec-
trum of C−1/ 2
G CLC−1/ 2
G , it is sensitive to local optima.
It is for instance rare that a dimension modelled by
a Gaussian is switched to the Parzen windows estima-
tor. Even though the algorithm will more easily switch
from the Parzen windows estimator to the Gaussian
model, it will typically stop too early, that is model
many dimensions using the Parzen windows estimator
rather than the better Gaussian. To solve these issues,
we propose an alternate algorithm, LCA-Gauss-Red,
which explores the space of dimensions modelled by a
Gaussian more aggressively using a search algorithm,
namely:
1. We run the algorithm LCA - Gauss for a few it-
erations (40 in our experiments);
2. We then “transfer” some columns from BL (the
Parzen windows model) to BG (the Gaussian
4
model), and rerun LCA - Gauss using these new
matrices as initialisations;
3. We iterate step 2 using a dichotomic search of the
optimal number of dimensions modelled by the
Gaussian, until a local optimum is found;
4. Once we have a locally optimum number of di-
mensions modelled by the Gaussian model, we run
LCA - Gauss to convergence.
5 Speeding up LCA
Computing the local covariance matrix of the points
using Eq. (3), (4) and (5) has a complexity in O(dn2 +
d2n + d3), with d the dimensionality of the data and
n the number of training points. Since this is imprac-
tical for large datasets, we can resort to sampling to
keep the cost linear in the number of datapoints. We
may further use low-rank or diagonal approximation to
achieve a complexity which grows quadratically with
d instead of cubically.
5.1 Averaging a subset of the local co-
variance matrices
Instead of averaging the local covariances over all dat-
apoints, we may only average them over a subset of
datapoints. This estimator is unbiased and, if the lo-
cal covariance matrices are not too dissimilar, which
is the assumption underlying LCA, then its variance
should remain small. This is equivalent to using a
minibatch procedure: every time we have a new mini-
batch of size B, we compute its local covariance ˆCL,
which is then averaged with the previously computed
CL using
CL ←γ
B
n CL + (1 −γ
B
n ) ˆCL (8)
to yield the updated CL. The exponent B/n is so
that γ, the discount factor, determines the weight of
the old covariance matrix after an entire pass through
the data, which makes it insensitive to the particular
choice of batch size. As opposed to many such al-
gorithms where the choice of γ is critical as it helps
retaining the information of previous batches, the lo-
cality of the EM estimate makes it less so. However, if
the number of datapoints used to estimate CL is not
much larger than the dimension of the data, we need
to set a higher γ to avoid degenerate covariance ma-
trices. In simulations, we found that using a value of
γ = . 6 worked well. Similarly, the size of the mini-
batch inﬂuences only marginally the ﬁnal result and
we found a value of 100 to be large enough.
5.2 Computing the local covariance
matrices using a subset of the dat-
apoints
Rather than using only a subset of local covariance ma-
trices, one may also wonder if using the entire dataset
to compute these matrices is necessary. Also, as the
number of datapoints grows, the chances of overﬁtting
increase. Thus, one may choose to use only a sub-
set of the datapoints to compute these matrices. This
will increase the local covariances, yielding a biased
estimate of the ﬁnal result, but may also act as a reg-
ulariser. In practice, for very large datasets, one will
want the largest neighbourhood size while keeping the
computational cost tractable.
Denoting ni the number of locations at which we esti-
mate the local covariance and nj the number of neigh-
bours used to estimate this covariance, the cost per
update is now O(d2[ni + nj] + dninj + d3). Since only
nj should grow with n, this is linear in the total num-
ber of datapoints.
Though they may appear similar, these are not “land-
mark” techniques (see, e.g., [11]) as there is still one
Gaussian component per datapoint, and the ni data-
points around which we compute the local covariances
are randomly sampled at every iteration.
6 Experiments
LCA has three main properties: ﬁrst, it transforms the
data to make it locally isotropic, thus being well-suited
for preprocessing the data before using a clustering
algorithm like spectral clustering; second, it extracts
relevant, non-Gaussian components in the data; third,
it provides us with a good density model through the
use of the Parzen windows estimator.
In the experiments, we will assess the performance of
the following algorithms: LCA, the original algorithm;
LCA-Gauss, using a multiplicative Gaussian compo-
nent, as described in Section 4; LCA-Gauss-Red, the
variant of LCA-Gauss using the more aggressive search
to ﬁnd a better number of dimensions to be modelled
by the Gaussian component. The MATLAB code for
LCA, LCA-Gauss and LCA-Gauss-Red is available at
http://nicolas.le-roux.name/code.html.
6.1 Improving clustering methods
We ﬁrst try to solve three clustering problems: one
for which the clusters are convex and the direc-
tion of interest does not have a Gaussian marginal
(Fig. (2), left), one for which the clusters are not con-
vex (Fig. (2), middle), and one for which the directions
5
of interest have almost Gaussian marginals (Fig. (2),
right). Following [1, 2], the data is progressively cor-
rupted by adding dimensions of white Gaussian noise,
then whitened. We compare here the clustering ac-
curacy, which is deﬁned as 100
n minP tr(EP ) where E
is the confusion matrix and P is the set of permuta-
tions over cluster labels, obtained with the following
ﬁve techniques:
1. Spectral clustering (SC) [19] on the whitened data
(using the code of [8]);
2. SC on the projection on the ﬁrst two components
found by FastICA using the best contrast function
and the correct number of components;
3. SC on the data transformed using the metric
learnt with LCA;
4. SC on the data transformed using the metric
learnt with the product of LCA and a Gaussian;
5. SC on the projection of the data found using
NGCA [4] with the correct number of compo-
nents.
Our choice of spectral clustering stems from its higher
clustering performance compared to K-means. Re-
sults are reported in Fig. (3). Because of the whiten-
ing, the Gaussian components in the ﬁrst dataset are
shrunk along the direction containing information. As
a result, even with little noise added, the information
gets swamped and spectral clustering fails completely.
On the other hand, LCA and its variants are much
more robust to the presence of irrelevant dimensions.
Though NGCA works very well on the ﬁrst dataset,
where there is only one relevant component, its per-
formance drops quickly when there are two relevant
components (note that, for all datasets, we provided
the true number of relevant dimensions as input to
NGCA). This is possibly due to the deﬂation proce-
dure which is not adapted when no single component
can be clearly identiﬁed in isolation. This is in contrast
with LCA and its variants which circumvent this issue,
thanks to their global optimisation procedure. Note
also that LCA-Gauss allows us to perform unsuper-
vised dimensionality reduction with the same perfor-
mance as previously proposed supervised algorithms
(e.g., [2]).
Figure (4) shows the clustering accuracy on the
three datasets for various numbers of EM iterations,
one iteration corresponding to Manifold Parzen Win-
dows [26] with a Gaussian kernel whose covariance ma-
trix is the data covariance kernel. As one can see, run-
ning the EM algorithm to convergence yields a signiﬁ-
cant improvement in clustering accuracy. The perfor-
mance of Manifold Parzen Windows could likely have
been improved with a careful initialisation of the orig-
inal kernel, but this would have been at the expense
of the simplicity of the algorithm.
6.2 LCA as a density model
We now assess the quality of LCA as a density model.
We build a density model of the USPS digits dataset,
a 256-dimensional dataset of handwritten digits. We
compared several algorithms:
•An isotropic Parzen windows estimator with the
bandwidth estimated using LCA (replacing Σ ∗ of
Eq. (5) by λI so that the two matrices have the
same trace);
•A Parzen windows estimator with diagonal metric
(equal to the diagonal of Σ ∗ in Eq. (5);
•A Parzen windows estimator with the full metric
as obtained using LCA;
•A single Gaussian model;
•A product of a Gaussian and a Parzen windows
estimator (as described in Section 4).
The models were trained on a set of 2000 datapoints
and regularized by penalizing the trace of Σ −1 (in the
case of the last model, both covariance matrices, local
and global, were penalized). The regularisation pa-
rameter was optimized on a validation set of 1000 dat-
apoints. For the last model, the regularisation param-
eter of the global covariance was set to the one yielding
the best performance for the full Gaussian model on
the validation set. Thus, we only had to optimize the
regularisation parameter for the local covariance.
The ﬁnal performance was then evaluated on a set of
3000 datapoints which had not been used for training
nor validation. We ran the experiment 20 times, ran-
domly selecting the training, validation and test set
each time.
Fig. (5) shows the mean and the standard error of the
negative log-likelihood on the test set. As one can see,
modelling all dimensions using the Parzen windows es-
timator leads to poor performance in high dimensions,
despite the regulariser and the leave-one-out criterion.
On the other hand, LCA-Gauss and LCA-Gauss-Red
clearly outperform all the other models, justifying our
choice of modelling some dimensions using a Gaussian.
Also, as opposed to the previous experiments, there is
no performance gain induced by the use of LCA-Gauss-
Red as opposed to LCA-Gauss, which we believe stems
from the fact that the switch from one model to the
other is easier to make when there are plenty of di-
mensions to choose from. The poor performance of
6
Figure 2: Noise-free data used to assess the robustness of K-means to noise. Left: mixture of two isotropic
Gaussians of unit variance and means [ −3, 0]⊤ and [3 , 0]⊤. Centre: two concentric circles with radii 1 and 2,
with added isotropic Gaussian noise of standard deviation . 1. Right: mixture of ﬁve Gaussians. The centre
cluster contains four times as many datapoints as the other ones.
0 20 40 6050
60
70
80
90
100
Noise dimensions
0 5 10 15 2050
60
70
80
90
100
Noise dimensions
0 20 40 60
30
40
50
60
70
80
90
100
Noise dimensions
Figure 3: Average clustering accuracy (100% = perfect clustering , chance is 50% for the ﬁrst two datasets and
20% for the last one) on 100 runs for varying number of dimensions o f noise added. The error bars represent one
standard error. Left: mixture of isotropic Gaussians presented in Fig. (2) (left). Centre: two concentric circles
presented in Fig. (2) (centre). Right: mixture of ﬁve Gaussians pr esented in Fig. (2) (right).
0 20 40 6050
60
70
80
90
100
Noise dimensions
 
 
Raw
MPW
LCA−10
LCA−50
LCA−200
0 5 10 15 2050
60
70
80
90
100
Noise dimensions
 
 
Raw
MPW
LCA−10
LCA−50
LCA−200
0 20 40 60
30
40
50
60
70
80
90
100
Noise dimensions
 
 
Raw
MPW
LCA−10
LCA−50
LCA−200
Figure 4: Average clustering accuracy (100% = perfect clustering , chance is 50% for the ﬁrst two datasets and
20% for the last one) on 100 runs for varying number of dimensions o f noise added and varying number of EM
iterations in the LCA algorithm (MPW = one iteration). The error bars represent one standard error. Left:
mixture of isotropic Gaussians presented in Fig. (2) (left). Centre : two concentric circles presented in Fig. (2)
(centre). Right: mixture of ﬁve Gaussians presented in Fig. (2) (r ight).
LCA-Full is a clear indication of the problems suﬀered
by Parzen windows in high dimensions.
6.3 Subsampling
We now evaluate the loss in performance incurred
by the use the subsampling procedure described in
Section 5, both on the train and test negative log-
7
N = 1000 N = 3000 N = 6000 N = 1000 N = 3000 N = 6000
B = 1000 6.43 ± 0.10 2.70 ± 0.06 − 0.10 ± 0.03 6.43 ± 0.10 2.80 ± 0.06 − 0.17 ± 0.02
B = 3000 6.58 ± 0.07 2.73 ± 0.05 − 0.06 ± 0.03 6.54 ± 0.07 2.80 ± 0.06 − 0.11 ± 0.02
B = 6000 6.22 ± 0.08 2.21 ± 0.03 0.00 ± 0.01 6.18 ± 0.07 1.98 ± 0.03 0.02 ± 0.02
N = 1000 N = 3000 N = 6000 N = 1000 N = 3000 N = 6000
B = 1000 2.12 ± 0.16 0.20 ± 0.08 0.65 ± 0.05 2.01 ± 0.17 0.07 ± 0.09 0.15 ± 0.03
B = 3000 2.11 ± 0.16 0.30 ± 0.07 0.46 ± 0.04 2.01 ± 0.17 0.16 ± 0.09 0.09 ± 0.03
B = 6000 1.54 ± 0.16 − 0.75 ± 0.07 − 0.01 ± 0.01 1.43 ± 0.15 − 1.07 ± 0.08 0.01 ± 0.02
Figure 6: Train (top) and test (bottom) negative log-likelihood diﬀer ences induced by the use of smaller batch
and neighbourhood sizes compared to the original model ( γ = 0, B = 6000, N = 6000) for γ = 0 . 3 (left) and
γ = 0 . 6 (right). A negative value means better performance.
LCA - Isotropic 269. 78 ±0. 18
LCA - Diagonal 109. 59 ±0. 56
LCA - Full 32. 98 ±0. 35
Gaussian 32. 27 ±0. 36
LCA - Gauss 19. 09 ±0. 39
LCA - Gauss - Red 19. 09 ±0. 39
Figure 5: Test negative log-likelihood on the USPS
digits dataset, averaged over 20 runs.
likelihoods. For that purpose, we used the USPS digit
recognition dataset, which contains 8298 datapoints in
dimension 256, which we randomly split into a train-
ing set of n = 6000 datapoints, using the rest as the
test set. We tested the following hyperparameters:
•Discount factor γ = 0 . 3, 0. 6, 0. 9 ,
•Batch size B = 1000 , 3000, 6000 ,
•Neighbourhood size N = 1000 , 3000, 6000 .
Fig. (6) show the log-likelihood diﬀerences induced
by the use of smaller batch sizes and neighbourhood
sizes. For each set of hyperparameters, 20 experiments
were run using diﬀerent training and test sets, and the
means and standard errors are reported. The results
for γ = 0 . 9 were very similar and are not included due
to space constraints.
Three observations may be made. First, reducing the
batchsize has little eﬀect, except when γ is small. Sec-
ond, reducing the neighbourhood size has a regular-
izing eﬀect at ﬁrst but drastically hurts the perfor-
mance if reduced too much. Third, the value of γ, the
discount factor, has little inﬂuence, but larger values
proved to yield more consistent test performance, at
the expense of slower convergence. The consistency of
these results shows that it is safe to use subsampling
(with values of γ = 0 . 6, B = 100 and N = 3000, for
instance) especially if the training set is very large.
7 Conclusion
Despite its importance, the learning of local or global
metrics is usually an overseen step in many practical
algorithms. We have proposed an extension of the gen-
eral bandwidth selection problem to the multidimen-
sional case, with a generalisation to the case where sev-
eral components are Gaussian. Additionally, we pro-
posed an approximate scheme suited to large datasets
which allows to ﬁnd a local optimum in linear time.
We believe LCA can be an important preprocessing
tool for algorithms relying on local distances, such as
manifold learning methods or many semi-supervised
algorithms. Another use would be to cast LCA within
the mean-shift algorithm, which ﬁnds the modes of the
Parzen windows estimator, in the context of image seg-
mentation [10]. In the future, we would like to extend
this model to the case where the metric is allowed to
vary with the position in space, to account for more
complex geometries in the dataset.
Acknowledgements
Nicolas Le Roux and Francis Bach are supported in
part by the European Research Council (SIERRA-
ERC-239993). We would also like to thank Warith
Harchaoui for its valuable input.
8
References
[1] F. Bach and Z. Harchaoui. Diﬀrac: a discrimi-
native and ﬂexible framework for clustering. In
Advances in Neural Information Processing Sys-
tems 20 , 2007.
[2] F. Bach and M. I. Jordan. Learning spectral clus-
tering. In Advances in Neural Information Pro-
cessing Systems 16 , 2004.
[3] Y. Bengio, H. Larochelle, and P. Vincent. Non-
local manifold Parzen windows. In Advances in
Neural Information Processing Systems 18 , pages
115–122. MIT Press, 2006.
[4] G. Blanchard, M. Kawanabe, M. Sugiyama,
V. Spokoiny, and K.-R. M¨ uller. In search of non-
Gaussian components of a high-dimensional dis-
tribution. J. Mach. Learn. Res. , 7:247–282, De-
cember 2006.
[5] J.M. Borwein and A.S. Lewis. Convex analysis
and nonlinear optimization, theory and examples,
2000.
[6] R. Boscolo, H. Pan, and V.P. Roychowdhury.
Independent component analysis based on non-
parametric density estimation. Neural Networks,
IEEE Transactions on , 15(1):55–65, 2004.
[7] A. Bowman. An alternative method of cross-
validation for the smoothing of density estimates.
Biometrika, 71(2):pp. 353–360, 1984.
[8] W.-Y. Chen, Y. Song, H. Bai, C.-J. Lin, and E. Y.
Chang. Parallel spectral clustering in distributed
systems. IEEE Transactions on Pattern Analysis
and Machine Intelligence , 33(3):568–586, 2011.
[9] Y.-S. Chow, S. Geman, and L.-D. Wu. Consistent
cross-validated density estimation. The Annals of
Statistics, 11(1):pp. 25–38, 1983.
[10] D. Comaniciu and P. Meer. Mean shift: A robust
approach toward feature space analysis. IEEE
Transactions on Pattern Analysis and Machine
Intelligence, 24:603–619, 2002.
[11] V. De Silva and J.B. Tenenbaum. Sparse multi-
dimensional scaling using landmark points. Tech-
nology, pages 1–41, 2004.
[12] A. P. Dempster, N. M. Laird, and D. B. Rubin.
Maximum likelihood from incomplete data via the
EM algorithm. Journal of the Royal Statistical
Society. Series B (Methodological) , 39(1):pp. 1–
38, 1977.
[13] R.P.W. Duin. On the choice of smoothing param-
eters for Parzen estimators of probability density
functions. Computers, IEEE Transactions on , C-
25(11):1175 –1179, 1976.
[14] T. Duong and M. Hazelton. Cross-validation
bandwidth matrices for multivariate kernel den-
sity estimation. Scandinavian Journal of Statis-
tics, 32(3):pp. 485–506, 2005.
[15] A. Globerson and S. Roweis. Metric learning
by collapsing classes. In Advances in Neural In-
formation Processing Systems , volume 18, pages
451–458. MIT Press, 2006.
[16] J. Goldberger, S. Roweis, G.E. Hinton, and
R. Salakhutdinov. Neighbourhood components
analysis. In Advances in Neural Information Pro-
cessing Systems 17 , 2004.
[17] A. Hyv¨ arinen, J. Karhunen, and E. Oja. Inde-
pendent component analysis . Wiley, New York,
2001.
[18] M. C. Jones and D. A. Henderson. Maximum
likelihood kernel density estimation: On the po-
tential of convolution sieves. Comput. Stat. Data
Anal., 53:3726–3733, August 2009.
[19] A. Ng, M. I. Jordan, and Y. Weiss. On spectral
clustering: Analysis and an algorithm. In Ad-
vances in Neural Information Processing Systems
14, 2001.
[20] B. Park and J. S. Marron. Comparison of data-
driven bandwidth selectors. Journal of the Amer-
ican Statistical Association , 85(409):pp. 66–72,
1990.
[21] E. Parzen. On estimation of a probability density
function and mode. The Annals of Mathematical
Statistics, 33(3):pp. 1065–1076, 1962.
[22] M. Rosenblatt. Remarks on some nonparametric
estimates of a density function. The Annals of
Mathematical Statistics , 27(3):pp. 832–837, 1956.
[23] M. Rudemo. Empirical choice of histograms and
kernel density estimators. Scandinavian Journal
of Statistics , 9(2):pp. 65–78, 1982.
[24] S. J. Sheather and M. C. Jones. A reliable data-
based bandwidth selection method for kernel den-
sity estimation. Journal of the Royal Statistical
Society. Series B (Methodological) , 53(3):pp. 683–
690, 1991.
[25] B. W. Silverman. Density estimation for statistics
and data analysis. 1986.
9
[26] P. Vincent and Y. Bengio. Manifold Parzen win-
dows. In Advances in Neural Information Pro-
cessing Systems 15 , pages 825–832. MIT Press,
2002.
[27] E. Xing, A. Ng, M. I. Jordan, and S. Russell. Dis-
tance metric learning with application to cluster-
ing with side-information. In Advances in Neural
Information Processing Systems 14 , 2002.
[28] L. Zelnik-Manor and P. Perona. Self-tuning spec-
tral clustering. In Advances in Neural Informa-
tion Processing Systems 17 , 2004.
Appendix
We prove here Proposition 1.
Proof If M1 is singular, then the minimum value is
−∞, because we can have B⊤
G M1BG bounded while
BGB⊤
G tends to + ∞(for example, if d1 = 1, and u1 is
such that M1u1 = 0, select BG = λu1 with λ →+∞).
The reasoning is similar for M2.
We thus assume that M1 and M2 are invert-
ible. We consider the eigendecomposition of
M−1/ 2
1 M2M−1/ 2
1 = U Diag(e)U⊤, which corresponds
to the generalized eigendecomposition of the pair
(M1, M 2).
Denoting A2 = U⊤M1/ 2
1 BL and A1 = U⊤M1/ 2
1 BG,
we have:
tr B⊤
G M1BG + tr B⊤
L M2BL −log det(BGB⊤
G + BLB⊤
L )
= tr A⊤
1 A1 + tr A⊤
2 Diag(e)A2
−log det(A1A⊤
1 + A2A⊤
2 ) + log det M1
= tr A⊤
1 A1 + tr A⊤
2 Diag(e)A2 −log det(A⊤
2 A2)
−log det(A⊤
1 (I −A2(A⊤
2 A2)−1A⊤
2 )A1) + log det M1 .
By taking derivatives with respect to A1, we get
A1 = ( I −Π 2)A1(A⊤
1 (I −Π 2)A1)−1 , (9)
with Π 2 = A2(A⊤
2 A2)−1A⊤
2 . By left-multiplying both
sides of Eq. (9) by A⊤
2 , we obtain
A⊤
2 A1 = 0 .
By left-multiplying by A⊤
1 , we get
A⊤
1 A1 = I .
Thus, we now need to minimize with respect to A2 the
following cost function
d1 + tr A⊤
2 Diag(e)A2 −log det(A⊤
2 A2) + log det M1
Let s be the vector of singular values of A2, ordered in
decreasing order and let the ei be ordered in increasing
order. We have:
tr Diag(e)A2A⊤
2 = −tr(−Diag(e)A2A⊤
2 ) ⩾
∑
i
eis2
i ,
with equality if and only if the eigenvectors of A2A⊤
2
are aligned with the ones of Diag( e) (the −ei being
also in decreasing order) (Theorem 1.2.1, [5]).
Thus, we have A2A⊤
2 = diag( s)2 with only d2 non-
zero elements in s. Let J2 be the index of non zero-
elements. We thus need to minimize
d1 + log det M1 +
∑
j∈J2
(ej s2
j −log s2
j ) ,
with optimum s2
j = e−1
j and value:
d1 + d2 + log det M1 +
∑
j∈J2
log ej .
Thus, we need to take J2 corresponding to the
smallest eigenvalues ej. If we also optimize with
respect to d2, then J2 must only contain the elements
smaller than 1.
10


Peer review:
[Review]:
First, we would like to thank the reviewers for their comments.

The main complaint was that the experiments were limited to toy problems. Since it is always hard to evaluate unsupervised learning algorithms (what is the metric of performance), the experiments were designed as a proof of concept. Hence, we agree with the reviewers and would love to see LCA tried and evaluated on real problems.

For the comment about the required modifications to avoid overfitting, there is truly only one parameter to set, i.e., the lambda parameter. All the others can easily be set to default values.