# An algorithmic Introduction to Clustering

Bernardo Gonzalez
Department of Computer Science and Engineering
UCSC
<EMAIL_ADDRESS>

The purpose of this document is to provide an easy introductory guide to
clustering algorithms. Basic knowledge and exposure to probability (random
variables, conditional probability, Bayes’ theorem, independence, Gaussian
distribution), matrix calculus (matrix and vector derivatives), linear algebra
and analysis of algorithm (specifically, time complexity analysis) is assumed.

Starting with Gaussian Mixture Models (GMM), this guide will visit different
algorithms like the well-known $k$-means, DBSCAN and Spectral Clustering (SC)
algorithms, in a connected and hopefully understandable way for the reader.
The first three sections (Introduction, GMM and $k$-means) are based on [1].
The fourth section (SC) is based on [2] and [5]. Fifth section (DBSCAN) is
based on [6]. Sixth section (Mean Shift) is, to the best of author’s
knowledge, original work. Seventh section is dedicated to conclusions and
future work.

Traditionally, these five algorithms are considered completely unrelated and
they are considered members of different families of clustering algorithms:

  * •

Model-based algorithms: the data are viewed as coming from a mixture of
probability distributions, each of which represents a different cluster. A
Gaussian Mixture Model is considered a member of this family

  * •

Centroid-based algorithms: any data point in a cluster is represented by the
central vector of that cluster, which need not be a part of the dataset taken.
$k$-means is considered a member of this family

  * •

Graph-based algorithms: the data is represented using a graph, and the
clustering procedure leverage Graph theory tools to create a partition of this
graph. Spectral Clustering is considered a member of this family

  * •

Density-based algorithms: this approach is capable of finding arbitrarily
shaped clusters, where clusters are defined as dense regions separated by low-
density regions. DBSCAN and Mean Shift are considered members of this family

In this document, we will try to present a more unified view of clustering, by
identifying the relationships between the algorithms mentioned. Some of the
results are not new, but they are presented in a cleaner, simpler and more
concise way. To the best of author’s knowledge, the interpretation of DBSCAN
as a climbing procedure, which introduces a theoretical connection between
DBSCAN and Mean shift, is a novel result.

Figure 1: What is the relation between these clustering algorihtms?

##  1 Introduction

The general version of the clustering problem is the following: given a set
$X$ of $n$ i.i.d. (independent and identically distributed) points
$x_{i}\in\mathcal{S},i=1,2,...,n$, and a distance function
$dist:\mathcal{S}\times\mathcal{S}\to\mathbb{R}_{+}$, cluster them in such a
way that points in the same group are less distant to each other than to those
in other clusters. In order to keep things simple, for the rest of this
document we will assume that $\mathcal{S}=\mathbb{R}^{d}$, but keep in mind
that most of the algorithms presented here can be generalized to spaces
different than $\mathbb{R}^{d}$. Also, a similarity measure
$sim:\mathcal{S}\times\mathcal{S}\to\mathbb{R}$ could be used instead of a
distance function. In this case, however, the goal is to cluster the points in
such a way that points in the same group are more similar to each other than
to those in other clusters.

###  1.1 Maximum Likelihood for the Gaussian distribution

Let $x\in\mathbb{R}^{d}$ be a $d$-dimensional random vector. A Gaussian
distribution (also known as Normal distribution) with mean vector
$\mu\in\mathbb{R}^{d}$ and covariance matrix $\Sigma\in\mathbb{R}^{d\times d}$
can be written in the form:

$\mathcal{N}(x|\mu,\Sigma)=\frac{1}{(2\pi)^{\frac{d}{2}}}\frac{1}{|\Sigma|^{\frac{1}{2}}}\exp{\Big{(}-\frac{1}{2}(x-\mu)^{T}\Sigma^{-1}(x-\mu)\Big{)}}$
(1)

where $|\Sigma|$ denotes the determinant of $\Sigma$ and $x^{T}$ is the
transpose of vector $x$.

Given a dataset $X$ composed by $n$ i.i.d. data points
$x_{1},x_{2},...,x_{n}$, where $\forall i,x_{i}\in\mathbb{R}^{d}$, we can use
a Gaussian Distribution as a model to fit this dataset. A common criterion for
fitting the parameters $\mu$ and $\Sigma$ of the Gaussian distribution given
an observed dataset is to find the values that maximize the likelihood
function:

$p(X|\mu,\Sigma)=\prod^{n}_{i=1}\mathcal{N}(x_{i}|\mu,\Sigma)=\prod^{n}_{i=1}\frac{1}{(2\pi)^{\frac{d}{2}}}\frac{1}{|\Sigma|^{\frac{1}{2}}}\exp\Big{(}-\frac{1}{2}(x_{i}-\mu)^{T}\Sigma^{-1}(x_{i}-\mu)\Big{)}$
(2)

This procedure is known as Maximum Likelihood Estimation (MLE). Given the form
of the Gaussian distribution, it is more convenient to maximize the $\log$ of
the likelihood function (maximizing this is equivalent because the logarithmic
function is a monotonically increasing function). Taking $\log$ on both sides

$\log
p(X|\mu,\Sigma)=-\frac{nd}{2}\log(2\pi)-\frac{n}{2}\log|\Sigma|-\frac{1}{2}\sum^{n}_{i=1}(x_{i}-\mu)^{T}\Sigma^{-1}(x_{i}-\mu)$
(3)

Maximizing with respect to $\mu$ (i.e., taking the derivative of $\log
p(X|\mu,\Sigma)$ w.r.t. $\mu$ and equalizing to zero), we obtain

$\mu_{MLE}=\frac{1}{n}\sum^{n}_{i=1}x_{i}$ (4)

Similarly, maximizing w.r.t. $\Sigma$

$\Sigma_{MLE}=\frac{1}{n}\sum^{n}_{i=1}(x_{i}-\mu_{MLE})(x_{i}-\mu_{MLE})^{T}$
(5)

As we can see, fitting one Gaussian distribution to a dataset $X$ using
maximum likelihood is a simple procedure, as shown in algorithm 1. The
complexity of this algorithm is $O(nd^{2})$, due the computation of matrix
$\Sigma_{MLE}$.

Input: Dataset $X=\\{x_{1},x_{2},...,x_{n}\\}$

Output: Learned values for $\mu$ and $\Sigma$

Estimate the parameters $\mu$ and $\Sigma$:

$\displaystyle\mu_{MLE}$
$\displaystyle\leftarrow\frac{1}{n}\sum^{n}_{i=1}x_{i}$
$\displaystyle\Sigma_{MLE}$
$\displaystyle\leftarrow\frac{1}{n}\sum^{n}_{i=1}(x_{i}-\mu_{MLE})(x_{i}-\mu_{MLE})^{T}$
return _$\mu_{MLE},\Sigma_{MLE}$_

Algorithm 1 Maximum Likelihood for one Gaussian Distribution

Figure 2(a) shows a dataset with 500 points in two dimensions. Figure 2(b)
shows the same dataset and the learned configuration for the Gaussian model.
The mean and the three standard-deviation contour for this model are shown as
a red point and a red ellipse, respectively.

(a) Dataset with 500 points in two dimensions

(b) Gaussian model learned from data using MLE

Figure 2: Fitting one Gaussian distribution to 500 data points. See text for
details

###  1.2 Sampling from the learned model

After fitting the model to the data distribution, we can sample from it and
compare how a set of samples looks with respect to the original dataset.
Figure 3 shows a set of 500 points sampled from the data (in red) along with
the 500 original data points (in blue). We can observe a couple of things: the
space covered by the sampled points is similar to the space covered by the
original data points, and any of the sampled points are "close" to at least
one of the original data points. Despite not being a quantitative measure,
figure 3 allow us to form an intuition of how well the model fitted the data
for this two dimensional dataset.

Figure 3: Points sampled from learned distribution (in red) along with
original dataset (in blue)

##  2 Gaussian Mixture Models (GMM)

###  2.1 Only one Gaussian: not flexible enough

(a) 1500 data points in two dimensions

(b) Gaussian model learned from data

(c) Points sampled from learned distribution (in red) along with original
dataset (in blue)

Figure 4: Fitting one Gaussian distribution to 1500 data points

As we can imagine, there exists datasets for which one Gaussian distribution
doesn’t offer the flexibility to fit the data points in the dataset to
analize. Figure 4(a) shows a dataset with 1500 points in two dimensions. We
will denote this dataset as $X_{1500}$. Figure 4(b) shows the dataset
$X_{1500}$ and the learned configuration for one Gaussian model. Like in the
previous example, the three standard-deviation contour is shown as a red
ellipse. However, in this case, we can observe how the model doesn’t capture
the data distribution as well as in the previous example. Specifically,
sampling from the learned distribution generates points which are "far" from
any of the original data points, as shown in figure 4(c). Also, we can observe
how the sampled points cover parts of the space not covered by the original
data points.

###  2.2 EM-GMM algorithm

One way to add more flexibility to our model is to use instead a Mixture of
Gaussians. A superposition of $k$ Gaussian densities of the form:

$p(x)=\sum^{k}_{c=1}\pi_{c}\mathcal{N}(x|\mu_{c},\Sigma_{c})$ (6)

is called a mixture of Gaussians. Each Gaussian density
$\mathcal{N}(x|\mu_{c},\Sigma_{c})$ is called a component of the mixture and
has its own mean $\mu_{c}$ and covariance matrix $\Sigma_{c}$. The parameters
$\pi_{c}\in\mathbb{R}$ are called mixing coefficients, and have two important
properties:

$\displaystyle 0\leq\pi_{c}\leq 1$ (7) $\displaystyle\sum^{k}_{c=1}\pi_{c}=1$
(8)

We can view the mixing coefficients $\pi_{c}$ as the prior probability $p(c)$
of picking the $c$-th component. Similarly, we can view the Gaussian densities
$\mathcal{N}(x|\mu_{c},\Sigma_{c})$ as the probability of $x$ conditioned on
$c$, $p(x|c)$ (i.e., the likelihood of point $x$ given the $c$-th component).
Using this notation, we can express equation (6) as:

$p(x)=\sum^{k}_{c=1}p(c)p(x|c)$ (9)

From here, we can see that we can compute the posterior probabilities $p(c|x)$
using Bayes’ theorem:

$\lambda_{i,c}\equiv
p(c|x_{i})=\frac{p(c)p(x_{i}|c)}{\sum^{k}_{j=1}p(j)p(x_{i}|j)}=\frac{\pi_{c}\mathcal{N}(x_{i}|\mu_{c},\Sigma_{c})}{\sum^{k}_{j=1}\pi_{j}\mathcal{N}(x_{i}|\mu_{j},\Sigma_{j})}$
(10)

Following the notation in [1], we denote the variables $\lambda_{i,c}$ as
"responsibilities", and we will use the following notation:
$\boldsymbol{\pi}\equiv\\{\pi_{1},...,\pi_{k}\\}$,
$\boldsymbol{\mu}\equiv\\{\mu_{1},...,\mu_{k}\\}$ and
$\boldsymbol{\Sigma}\equiv\\{\Sigma_{1},...,\Sigma_{k}\\}$. These
responsibilities $\lambda_{i,c}$ can be interpreted as a measure of how much
the point $x_{i}$ "belongs" to the $c$-th component, given the current
estimates $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$. Similar to what we did
with one Gaussian distribution, we can try to fit the parameters of this model
given an observed dataset $X$ using maximum likelihood. Using the i.i.d.
assumption and equation (6) we have

$\log
p(X|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})=\log\Big{(}\prod^{n}_{i=1}p(x_{i}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})\Big{)}=\sum^{n}_{i=1}\log\Big{(}\sum^{k}_{c=1}\pi_{c}\mathcal{N}(x_{i}|\mu_{c},\Sigma_{c})\Big{)}$
(11)

The situation now is much more complex than with a single Gaussian, due to the
presence of the summation inside the logarithm. As a result, the maximum
likelihood solution for the parameters of the GMM no longer has a closed-form
analytical solution. Maximizing equation (11) with respect to $\mu_{c}$, we
obtain

$\mu_{c}=\frac{1}{n_{c}}\sum^{n}_{i=1}\lambda_{i,c}x_{i}$ (12)

where $n_{c}=\sum^{n}_{i=1}\lambda_{i,c}$. Similarly, maximizing equation (11)
w.r.t. $\Sigma_{c}$

$\Sigma_{c}=\frac{1}{n_{c}}\sum^{n}_{i=1}\lambda_{i,c}(x_{i}-\mu_{c})(x_{i}-\mu_{c})^{T}$
(13)

Finally, we need to maximize equation (11) w.r.t. $\pi_{c}$. However, in this
case we need to take into acount the constraints (7) and (8). This can be
achieved using Lagrange multipliers, getting as result

$\pi_{c}=\frac{n_{c}}{n}$ (14)

Despite looking as closed-form solutions, equations (12),(13) and (14) are not
such, because all of them are functions of the responsibilities
$\lambda_{i,c}$. However, these responsibilities depends on
$\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}$ in a complex way, as
we can see from equation (10).

However, the last three equations suggest a simple iterative algorithm to find
a solution to the GMM maximum likelihood problem:

  * •

First, we choose some initial values for the means $\mu_{c}$, matrices
$\Sigma_{c}$ and mixing coefficients $\pi_{c}$

  * •

Then, we alternate between the following two updates, that we shall call the E
step (from Expectation) and M step (from Maximization):

    * –

In the E step, we use the current values of $\mu_{c}$, $\Sigma_{c}$ and
$\pi_{c}$ to compute the responsibilities $\lambda_{i,c}$ using equation (10)

    * –

Then, we use these $\lambda_{i,c}$ in the M step to re-estimate the values of
$\mu_{c}$, $\Sigma_{c}$ and $\pi_{c}$ using equations (12),(13) and (14),
respectively

It is important to highlight that when computing the update for
$\Sigma^{new}_{c}$, we use the updated values of $\mu^{new}_{c}$. The
algorithm iterates until a convergence criteria is met. In practice, the
algorithm is deemed to have converged when the change in the log likelihood
function , or alternatively in the parameters, falls below some threshold.

This iterative procedure using the E step and M step is an instance of a more
general algorithm called Expectation-Maximization algorithm. The EM algorithm
is an iterative procedure to find maximum likelihood solutions for models
having latent variables. It should be emphasized that commonly, the log
likelihood function has multiple local maxima, and that the EM algorithm is
not guaranteed to find the global maxima. The EM algorithm for Gaussian
Mixtures models (EM-GMM algorithm) is shown in algorithm 2.

EM-GMM is the first clustering algorithm presented in this document: it takes
a collection $X$ of $n$ points in a $d$ dimensional space and group them into
$k$ different clusters $\Gamma_{1},\Gamma_{2},...,\Gamma_{k}$ using a mixture
of Gaussians as a model. Every Gaussian distribution models a cluster
$\Gamma_{c}$. Figure 5(a) shows how a GMM model with three components is able
to fit the previous dataset more tightly than the one-Gaussian model. Also,
sampling from the GMM model generates points which are more similar to points
in the original dataset, as shown in figure 5(b). We can also observe how the
space covered by the sampled points is now more similar to the space covered
by the original data points.

(a) GMM with three components learned from data

(b) Points sampled from GMM (in red) along with original dataset (in blue)

Figure 5: Fitting a GMM with three components to 1500 data points

Input: Dataset $X$, number of components $k$, and initial values for every
$\pi_{c},\mu_{c},\Sigma_{c}$

Output: Learned values for every $\pi_{c},\mu_{c},\Sigma_{c}$ and
$\lambda_{i,c}$

while  _Convergence criterion is not satisfied_ do

E step. Evaluate the responsibilities $\lambda_{i,c}$ using the current values
of $\pi_{c},\mu_{c},\Sigma_{c}$:

$\lambda_{i,c}\leftarrow\frac{\pi_{c}\mathcal{N}(x_{i}|\mu_{c},\Sigma_{c})}{\sum^{k}_{j=1}\pi_{j}\mathcal{N}(x_{i}|\mu_{j},\Sigma_{j})}\;$

M step. Re-estimate the parameters using the current responsibilities
$\lambda_{i,c}$:

$\displaystyle n_{c}$ $\displaystyle\leftarrow\sum^{n}_{i=1}\lambda_{i,c}$
$\displaystyle\mu^{new}_{c}$
$\displaystyle\leftarrow\frac{1}{n_{c}}\sum^{n}_{i=1}\lambda_{i,c}x_{i}$
$\displaystyle\Sigma^{new}_{c}$
$\displaystyle\leftarrow\frac{1}{n_{c}}\sum^{n}_{i=1}\lambda_{i,c}(x_{i}-\mu^{new}_{c})(x_{i}-\mu^{new}_{c})^{T}$
$\displaystyle\pi^{new}_{c}$ $\displaystyle\leftarrow\frac{n_{c}}{n}\;$

return
_$\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma},\boldsymbol{\lambda}$_

Algorithm 2 EM-GMM

###  2.3 EM-GMM complexity

At every iteration, the two updates more computationally intensive are the
update of $\lambda_{i,c}$ in E step and the update of the matrices
$\Sigma_{c}$ in the M step.

Specifically, at every iteration, the E step requires the evaluation of $nk$
Gaussian distributions, which involves the computation (for every cluster
$\Gamma_{c}$) of the determinant and the inverse of a matrix of size $d\times
d$, both with complexity $O(d^{2.807})$ for practical purposes (using Strassen
algorithm for matrix multiplication, however there are algorithms
asymptotically faster). Also, for every $x_{i}$, we need to compute a
quadratic form w.r.t. every cluster $\Gamma_{c}$, with complexity $O(d^{2})$,
therefore resulting in an overall complexity of $O(nd^{2}k)$. Therefore, the
overall complexity for E step at every iteration is $O(nd^{2}k+d^{2.807}k)$.

Similarly, the update of any of the matrices $\Sigma_{c}$ requires the
computation and addition of $n$ matrices of size $d\times d$, with complexity
$O(nd^{2})$, therefore resulting in an overall complexity of $O(nd^{2}k)$ for
M step at every iteration.

Therefore, the complexity of the EM-GMM algorithm is
$O(nd^{2}kT+d^{2.807}kT)$, where $T$ denotes the number of iterations. This
complexity and other important properties of the EM-GMM algorithm are
summarized in table 1.

Table 1: EM-GMM algorithm | GMM
---|---
Number of hyperparameters | 1
Hyperparameters | number of components $k$
Time complexity | $O(nd^{2}kT+d^{2.807}kT)$
Outlier detection | No
Type of assignment | soft
distance/similarity measure | Mahalanobis distance (See $k$-means section for details)

###  2.4 GMM clustering

A GMM returns what is called a soft assignment for every data point $x_{i}$:
instead of assigning it to a specific cluster $\Gamma_{c}$, it returns a
"responsibility" $\lambda_{i,c}$, expressing how "responsible" is cluster
$\Gamma_{c}$ for data point $x_{i}$. The set of $\lambda_{i,c}$ for a specific
point $x_{i}$ forms a distribution over the clusters
$\Gamma_{1},\Gamma_{2},...,\Gamma_{k}$. On the other hand, a hard assignment
is an assignment where every point $x_{i}$ is assigned to exactly one cluster
$\Gamma_{c}$.

One simple way to turn the soft assignment returned by a GMM into a hard
assignment is to assign point $x_{i}$ to the cluster $\Gamma_{c}$ with larger
responsability $\lambda_{i,c}$. Using this rounding method, we can observe in
figure 6 how the clusters found by the GMM model (right side) compare to the
labels from the $X_{1500}$ dataset (left side).

Figure 6: Clustering induced by ground truth (labels) and clustering learnt by
GMM

####  2.4.1 Clustering quality measures: Adjusted Mutual Information (AMI) and
Adjusted Rand Index (ARI)

In order to present a more quantitative measure of how similar are the labels
from the data and the clusters found by the GMM, we introduce two widely-used
methods for comparing a pair of clusterings (or more generally, two partitions
of a set): the Adjusted Mutual Information (AMI) [11] and the Adjusted Rand
Index (ARI) [12]. AMI is an entropy-based measure that quantifies the overlap
between the clusters found and the clusters defined by the labels. Similarly,
ARI is a measure that quantifies the number of agreements and disagreements
between the clusters found and the clusters defined by the labels. Both scores
take the value of 1 when two clusterings are identical. For more details,
please refer to [11] and [12], respectively.

The AMI and the ARI scores between the labels and the clusters found by the
GMM for the dataset shown in figure 6 are 0.992 and 0.996, respectively.

##  3 $k$-means

$k$-means is another clustering algorithm, with some similarities with the EM-
GMM algorithm (both are iterative algorithms and both need the number of
clusters $k$ as an input). Can we derive the $k$-means algorithm from EM-GMM?
In the next section we’ll observe how we can derive the $k$-means algorithm as
a particular limit of the EM-GMM algorithm.

###  3.1 From GMM to $k$-means

Consider a GMM model with all covariance matrices equal to the identity matrix
$I$ times a variance parameter $\epsilon$, i.e., $\Sigma_{c}=\epsilon I$
$\forall c$. Then, for point $x$ and cluster $\Gamma_{c}$

$p(x|\mu_{c},\Sigma_{c})=\frac{1}{(2\pi\epsilon)^{\frac{d}{2}}}\exp{\Big{(}-\frac{1}{2\epsilon}\left\lVert
x-\mu_{c}\right\rVert^{2}\Big{)}}$ (15)

where we used the identity $\left\lVert x\right\rVert^{2}=x^{T}x$. Lets now
consider algorithm 2 for a GMM with $k$ Gaussians of this form and treat
$\epsilon$ as a fixed constant. In the E step, the responsibilities
$\lambda_{i,c}$ for data point $x_{i}$ are now given by

$\lambda_{i,c}=\frac{\pi_{c}\exp{\Big{(}-\frac{1}{2\epsilon}\left\lVert
x_{i}-\mu_{c}\right\rVert^{2}\Big{)}}}{\sum^{k}_{j=1}\pi_{j}\exp{\Big{(}-\frac{1}{2\epsilon}\left\lVert
x_{i}-\mu_{j}\right\rVert^{2}\Big{)}}}$ (16)

If we consider the limit $\epsilon\to 0$, we observe that in the denominator
the term for which $\left\lVert x_{i}-\mu_{j}\right\rVert^{2}$ is smallest
will go to zero most slowly, and hence the responsibilities $\lambda_{i,c}$
for the data point $x_{i}$ all go to zero except for term $j$, for which the
responsibility $\lambda_{i,j}$ will go to one. Note that this holds
independently of the values of $\pi_{c}$ as long as none of the $\pi_{c}$ is
zero. Therefore, in the limit $\epsilon\to 0$, algorithm 2 returns a hard
assignment: it assigns every data point $x_{i}$ to just one specific cluster
$\Gamma_{c}$. It is easy to observe that data points will be assigned to the
cluster having the closest mean $\mu_{c}$. This is the E step in the $k$-means
algorithm: assign every point to the cluster with the current closest mean
$\mu_{c}$. Also, instead of using $\lambda_{i,c}$ as the assignments notation,
we will follow [1] and use $r_{i,c}$ instead as assignments notation in the
$k$-means algorithm.

Continuing with the EM-GMM algorithm, lets now observe what changes in the M
step. The first thing that we can observe is that the last two equations
(updating the covariance matrices $\Sigma^{new}_{c}$ and $\pi^{new}_{c}$) are
no longer needed: all matrices $\Sigma_{c}$ are constant and equal, and we no
longer need $\pi_{c}$ for the E step. Also, we can merge the first two
equations in M step into one equation which updates the means $\mu_{c}$ to the
mean of the current cluster (after the E step update):

$\mu_{c}=\frac{\sum^{n}_{i=1}r_{i,c}x_{i}}{\sum^{n}_{i=1}r_{i,c}}$ (17)

So, the answer to the question opening this section "Can we derive the
$k$-means algorithm from the EM-GMM one?" is yes: using the EM-GMM algorithm
with a GMM with all covariance matrices $\Sigma_{c}=\epsilon I$ and taking the
limit when $\epsilon\to 0$, leads to the $k$-means algorithm, shown in
algorithm 3. The consequences of these changes made to the EM-GMM algorithm is
that the $k$-means algorithm lose some flexibility w.r.t. a GMM: all points
are assigned now to a single cluster, and, because all covariances matrices
are diagonal matrices, all data points are assigned to the cluster with the
closest mean (i.e., for $k$-means, "closest" means closest w.r.t. to Euclidean
distance, while for a GMM, "closest" means closest w.r.t. to Mahalanobis
distance [13]).

However, the $k$-means algorithm is more efficient than EM-GMM: we do not need
to evaluate the probabilities $\mathcal{N}(x_{i}|\mu_{c},\Sigma_{c})$ for
every data point and every cluster at E step. Also, is not necessary to update
the covariance matrices $\Sigma^{new}_{c}$ at M step. Avoiding these updates
reduces significantly the complexity of the algorithm, which is $O(ndkT)$,
because at E step we need to compute the distance between every data point
$x_{i}$ and every cluster mean $\mu_{c}$, at every iteration. This complexity
and other important properties of the $k$-means algorithm are summarized in
table 2.

Figure 7 illustrates the main difference between $k$-means and a GMM. We can
observe how $k$-means is not flexible enough to capture the elliptical shape
of the data, in contrast to the results obtained using a GMM. The main reason
for this is the distance function used by each algorithm, the Mahalanobis
distance for GMMs and the Euclidean distance for $k$-means. The circles shown
for $k$-means are for illustration purpose only: as mentioned before,
$k$-means does not compute the covariance matrices $\Sigma_{c}$.

(a) A GMM is flexible enough to fit the data

(b) $k$-means lacks the flexibility to fit this dataset

Figure 7: Difference between a GMM and $k$-means

Input: Dataset $X$, number of clusters $k$, and initial values for every
$\mu_{c}$

Output: Learned values for every $\mu_{c}$ and $r_{i,c}$

while  _Convergence criterion is not satisfied_ do

E step. Using the current values of $\mu_{c}$, assign every point to its
closest cluster:

$r_{i,c}\leftarrow\begin{cases}1&\text{if
$c=argmin_{j}||x_{i}-\mu_{j}||^{2}$}\\\ 0&\text{otherwise}\end{cases}$

M step. Re-estimate the cluster means using the current assignations
$r_{i,c}$:

$\mu_{c}\leftarrow\frac{\sum^{n}_{i=1}r_{i,c}x_{i}}{\sum^{n}_{i=1}r_{i,c}}$

return _$\boldsymbol{\mu},\boldsymbol{r}$_

Algorithm 3 $k$-means Table 2: $k$-means algorithm |  $k$-means
---|---
Number of hyperparameters | 1
Hyperparameters | number of clusters $k$
Time complexity | $O(ndkT)$
Outlier detection | No
Type of assignment | hard
distance/similarity measure | Euclidean distance

###  3.2 $k$-means as an optimization problem

Before continuing, we’ll introduce a different formulation for the $k$-means
problem that will be useful in the next section. Given a dataset $X$ with $n$
data points $x_{1},x_{2},...,x_{n}$, the $k$-means algorithm tries to find $k$
points in the space (the means $\mu_{c}$) and assign every data point $x_{i}$
to one of these $k$ points such that the distance between $x_{i}$ and
$\mu_{c}$ is minimized. We can then think of the $k$-means algorithm as an
iterative algorithm that returns a local minima for the following optimization
problem:

$\displaystyle\min_{r_{i,c},\mu_{c}}$
$\displaystyle\sum^{k}_{c=1}\sum_{i\in\Gamma_{c}}\left\lVert
x_{i}-\mu_{c}\right\rVert^{2}$ (18) s.t. $\displaystyle
r_{i,c}\in\\{0,1\\}\qquad(\forall i=1,2,...n;c=1,2,...k)$

where we have used the notation $i\in\Gamma_{c}$ to denote that point $x_{i}$
belongs to cluster $\Gamma_{c}$, i.e., $r_{i,c}=1$ and $r_{i,j}=0$ $\forall
j\neq c$. We’ll denote this optimization problem as the vanilla $k$-means
problem. A solution for this problem consists of $k$ cluster means $\mu_{c}$
and a hard assignment $r_{i,c}$ for every data point $x_{i}$. The $k$-means
algorithm returns a local minima for this optimization problem using the EM
algorithm.

##  4 Spectral Clustering (SC)

Spectral Clustering is a family of algorithms that use the eigenvectors of a
matrix derived from a dataset $X$ to cluster the data points in $X$. Each
specific algorithm use the eigenvectors in a slightly different way. In this
document, we will focus on the spectral algorithm of Ng, Jordan and Weiss [3]
that we’ll denote as the NJW-SC algorithm.

###  4.1 An extension to $k$-means: Weighted kernel $k$-means

A major drawback of the $k$-means algorithm is that it can not separate
clusters that are non-linearly separable in input space $\mathcal{S}$ (i.e.,
if the intersection of the convex hull of the clusters is not empty). One
approach for tackling this problem is using the kernel method, leading to an
algorithm called kernel $k$-means.

####  4.1.1 The kernel method

The kernel method is commonly used in different Machine Learning algorithms.
The basic idea behind the kernel method is to map each data point into a
different space $\mathcal{F}$ (usually a high dimensional space) via a map
$\phi:\mathcal{S}\to\mathcal{F}$, and then use the learning algorithm (in our
case, the $k$-means algorithm) in this space $\mathcal{F}$. In our specific
case, we would like to have a map $\phi$ such that in the space $\mathcal{F}$
the data points are now linearly separable. Figure 8 shows a simple example of
the kernel method. The concentric circles dataset shown in figure 8(a) is not
linearly separable in the original (two dimensional) space $\mathcal{S}$, but
using the mapping $\phi(x_{i})=[x_{i},\left\lVert x_{i}\right\rVert]$ (i.e.,
using the distance from any point $x_{i}$ to the origin as its third
coordinate in the three dimensional space $\mathcal{F}$), now the circles are
linearly separable in $\mathcal{F}$, as shown in figures 8(b) and 8(c).

(a) Concentric circles dataset

(b) Concentric circles dataset in space $\mathcal{F}$

(c) Hyperplane separating the concentric circles dataset in $\mathcal{F}$

Figure 8: An example of the kernel method

One of the key ideas behind the kernel method is to realize that most distance
functions (like the Euclidean distance) can be expressed in terms of dot
products. Therefore, even if now the data points lies in a possibly high
dimensional space, if we can compute the dot product between any pair of
points in an efficient manner, then we can compute the distance between them
also in an efficient way. This is exactly what a kernel function is: a
function $\kappa:\mathcal{S}\times\mathcal{S}\to\mathbb{R}$ which computes the
dot product between a couple of (mapped) points $\phi(x_{i})$ and
$\phi(x_{j})$ without explicitly computing the mappings $\phi(x_{i})$ and
$\phi(x_{j})$, i.e., $\kappa(x_{i},x_{j})=\phi^{T}(x_{i})\phi(x_{j})$. Popular
kernel functions are the polynomial kernel
$\kappa_{P}(x_{i},x_{j})=(x^{T}_{i}x_{j}+c)^{b}$, $c,b\in\mathbb{R}$ or the
Gaussian kernel $\kappa_{G}(x_{i},x_{j})=\exp{\Big{(}-\frac{\left\lVert
x_{i}-x_{j}\right\rVert^{2}}{2\sigma^{2}}\Big{)}}$, $\sigma\in\mathbb{R}$.

####  4.1.2 Weighted kernel $k$-means

Another extension to the $k$-means algorithm is to associate each data point
$x_{i}$ with a weight $w_{i}\in\mathbb{R}$. Depending on the context, weight
$w_{i}$ usually express how important point $x_{i}$ is. Adding both the kernel
method and the weights to the optimization formulation of the $k$-means
problem (18) results in the following optimization problem:

$\displaystyle\min_{r_{i,c}}$
$\displaystyle\sum^{k}_{c=1}\sum_{i\in\Gamma_{c}}w_{i}\left\lVert\phi(x_{i})-m_{c}\right\rVert^{2}$
(19) s.t. $\displaystyle r_{i,c}\in\\{0,1\\}\qquad(\forall
i=1,2,...n;c=1,2,...k)$

Note that weights $w_{i}$ are not optimization variables, they are known
constants, an input to the optimization problem. Also, observe that we have
used $m_{c}$ to denote the weighted mean of the points in cluster
$\Gamma_{c}$, but in the space $\mathcal{F}$:

$\displaystyle m_{c}=\frac{\sum_{i\in\Gamma_{c}}w_{i}\phi(x_{i})}{s_{c}}$ (20)

where $s_{c}=\sum_{i\in\Gamma_{c}}w_{i}$ is the total weight of cluster
$\Gamma_{c}$. Observe that $m_{c}$ is the point that minimizes

$\sum_{i\in\Gamma_{c}}w_{i}\left\lVert\phi(x_{i})-z\right\rVert^{2}$

w.r.t. $z$, i.e., $m_{c}$ is the "best" cluster representative for cluster
$\Gamma_{c}$ in space $\mathcal{F}$. We can rewrite the objective function of
problem (19) as follows

$\displaystyle\sum^{k}_{c=1}\sum_{i\in\Gamma_{c}}w_{i}\left\lVert\phi(x_{i})-m_{c}\right\rVert^{2}$
$\displaystyle=\sum^{k}_{c=1}\sum_{i\in\Gamma_{c}}w_{i}\big{(}\phi(x_{i})-m_{c}\big{)}^{T}\big{(}\phi(x_{i})-m_{c}\big{)}$
$\displaystyle=\sum^{k}_{c=1}\sum_{i\in\Gamma_{c}}w_{i}\big{(}\phi^{T}(x_{i})\phi(x_{i})-2\phi^{T}(x_{i})m_{c}+m^{T}_{c}m_{c}\big{)}$
$\displaystyle=\sum^{k}_{c=1}\Big{(}\sum_{i\in\Gamma_{c}}w_{i}\phi^{T}(x_{i})\phi(x_{i})-2\sum_{i\in\Gamma_{c}}w_{i}\phi^{T}(x_{i})m_{c}+\sum_{i\in\Gamma_{c}}w_{i}m^{T}_{c}m_{c}\Big{)}$
$\displaystyle=\sum^{k}_{c=1}\Big{(}\sum_{i\in\Gamma_{c}}w_{i}\phi^{T}(x_{i})\phi(x_{i})-2s_{c}m^{T}_{c}m_{c}+s_{c}m^{T}_{c}m_{c}\Big{)}$
$\displaystyle=\sum^{k}_{c=1}\Big{(}\sum_{i\in\Gamma_{c}}w_{i}\phi^{T}(x_{i})\phi(x_{i})-s_{c}m^{T}_{c}m_{c}\Big{)}$
$\displaystyle=\sum^{k}_{c=1}\sum_{i\in\Gamma_{c}}w_{i}\phi^{T}(x_{i})\phi(x_{i})-\sum^{k}_{c=1}s_{c}m^{T}_{c}m_{c}$

where fourth equality is due equation (20). Also, note that the first term in
the last equation does not depend on the assignments $r_{i,c}$ or the cluster
centers $m_{c}$, so we can remove it from the objective function, obtaining
the following objective function

$-\sum^{k}_{c=1}s_{c}m^{T}_{c}m_{c}=-\sum^{k}_{c=1}\frac{\sum_{i,j\in\Gamma_{c}}w_{i}w_{j}\phi^{T}(x_{i})\phi(x_{j})}{s_{c}}$
(21)

From the form of this equation, we can see that the quantity to minimize
(across all the clusters) is the negative sum of the weighted dot product (in
space $\mathcal{F}$) between all pair of points in the same cluster,
normalized by the total weight of the cluster. In other words, and
interpreting the dot product as a similarity measure (in terms of the
direction of the data points), we are trying to maximize the normalized
similarity among points in the same cluster, where the normalization term is
the weight of the cluster, and the similarity function between two points is
defined by the kernel $\kappa()$.

Letting $W\in\mathbb{R}^{n\times n}$ be the diagonal matrix with all the
weights $w_{i}$ in the diagonal,
$W_{c}\in\mathbb{R}^{|\Gamma_{c}|\times|\Gamma_{c}|}$ (where $|\Gamma_{c}|$
denotes the number of elements in cluster $\Gamma_{c}$) be the diagonal matrix
of the weights in cluster $\Gamma_{c}$,
$\Phi\in\mathbb{R}^{dim(\mathcal{F})\times n}$ (where $dim(\mathcal{F})$
denotes the dimension of space $\mathcal{F}$) be the matrix formed by
horizontally concatenating the points $\phi(x_{i})$ (i.e.,
$\Phi=[\phi(x_{1}),\phi(x_{2}),...,\phi(x_{n})]$),
$\Phi_{c}\in\mathbb{R}^{dim(\mathcal{F})\times|\Gamma_{c}|}$ the matrix formed
by horizontally concatenating the points $\phi(x_{i})$ belonging to cluster
$\Gamma_{c}$ and $e_{c}$ the vector of ones of size $|\Gamma_{c}|$, we can
rewrite $m_{c}$ as

$\displaystyle m_{c}=\frac{\Phi_{c}W_{c}e_{c}}{s_{c}}$ (22)

Using this equation, we can rewrite equation (21) as

$\displaystyle-\sum^{k}_{c=1}s_{c}m^{T}_{c}m_{c}$
$\displaystyle=-\sum^{k}_{c=1}\frac{e^{T}_{c}W_{c}\Phi^{T}_{c}\Phi_{c}W_{c}e_{c}}{s_{c}}$
$\displaystyle=-\sum^{k}_{c=1}\frac{e^{T}_{c}W_{c}\Phi^{T}_{c}}{\sqrt{s_{c}}}\frac{\Phi_{c}W_{c}e_{c}}{\sqrt{s_{c}}}$
$\displaystyle=-Tr(Y^{T}W^{\frac{1}{2}}\Phi^{T}\Phi W^{\frac{1}{2}}Y)$
$\displaystyle=-Tr(Y^{T}W^{\frac{1}{2}}KW^{\frac{1}{2}}Y)$

where $Tr(A)$ denotes the trace of matrix $A$, $K=\Phi^{T}\Phi$ is the kernel
matrix of the data, i.e., $K_{ij}=\kappa(x_{i},x_{j})$, and

$Y=\begin{bmatrix}\frac{W^{\frac{1}{2}}_{1}e_{1}}{\sqrt{s_{1}}}&0&0&0\\\
0&\frac{W^{\frac{1}{2}}_{2}e_{2}}{\sqrt{s_{2}}}&0&0\\\ 0&0&...&0\\\
0&0&0&\frac{W^{\frac{1}{2}}_{k}e_{k}}{\sqrt{s_{k}}}\\\ \end{bmatrix}$ (23)

is the assignment matrix, where all the zeros denotes zero matrices of
appropiate size. Note that $Y$ is an $n\times k$ orthonormal matrix, i.e.,
$Y^{T}Y=I$. We can then re-cast the optimization problem (19) as

$\displaystyle\max_{Y\in\nu_{Y}^{n\times k}}$ $\displaystyle
Tr(Y^{T}W^{\frac{1}{2}}KW^{\frac{1}{2}}Y)$ (24) s.t. $\displaystyle Y^{T}Y=I$

where
$\nu_{Y}=\\{0\\}\bigcup\left\\{\sqrt{\frac{w_{i}}{s_{c}}}\right\\},i\in\Gamma_{c}$.
We will denote problem (24) as the matrix version of the weighted kernel
$k$-means problem. There are two important things to observe from problem
(24). The first one is that, if we set all the weights $w_{i}$ equal to one
and set the function $\phi()$ to be the identity function, then problem (19)
turns into problem (18), the vanilla $k$-means problem, and the corresponding
matrix version of this problem turns into the following optimization problem

$\displaystyle\max_{Y\in\nu_{Y}^{n\times k}}$ $\displaystyle Tr(Y^{T}GY)$ (25)
s.t. $\displaystyle Y^{T}Y=I$

where
$\nu_{Y}=\\{0\\}\bigcup\left\\{\sqrt{\frac{1}{s_{c}}}\right\\},i\in\Gamma_{c}$,
and $G$ is the Gram matrix of the data points, i.e., $G_{ij}=x^{T}_{i}x_{j}$.
We will denote problem (25) as the matrix version of the vanilla $k$-means
problem.

The second thing is that, just like the $k$-means algorithm can find a local
minima for the vanilla $k$-means problem, we can use a modified version of the
$k$-means algorithm to find a local minima for the weighted kernel $k$-means
problem. This is shown in algorithm 4. There is something important to
highlight with this algorithm. Following it naively will get you into a
situation where you need the mapped points $\phi(x_{i})$ in both steps E and
M. However, as we mentioned previously, the whole idea of the kernel method is
to avoid the explicit computation of $\phi(x_{i})$, and instead use only the
dot products $\phi^{T}(x_{i})\phi(x_{j})=\kappa(x_{i},x_{j})$. So, for step E,
when computing the squared distance $||\phi(x_{i})-m_{j}||^{2}$, it is implied
that we need to use the representation of the mean $m_{j}$ using equation
(20), substitute it in the equation of the squared distance, re-write this
equation as a dot product, and then use the distributive property of the dot
product to express $||\phi(x_{i})-m_{j}||^{2}$ as a sum of weighed dot
products $\phi^{T}(x_{i})\phi(x_{k})$ that we can compute using the kernel
function $\kappa()$. Finally, observe that we don’t need to compute the value
of $m_{c}$ in the M step, because, as mentioned previously, it can be only
implicitly computed via $\kappa()$ (we can’t explicitly compute it, because we
don’t have access to the mapped points $\phi(x_{i})$). This is the reason why,
in contrast to the vanilla $k$-means algorithm, the weighted kernel $k$-means
algorithm only returns the assignments $r_{i,c}$.

Input: Dataset $X$, number of clusters $k$, weights $w_{i}$, kernel function
$\kappa()$ and initial values for every $m_{c}$

Output: Learned values for $r_{i,c}$

while  _Convergence criterion is not satisfied_ do

E step. Using the current values of $m_{c}$, assign every point to its closest
cluster:

$r_{i,c}\leftarrow\begin{cases}1&\text{if
$c=argmin_{j}||\phi(x_{i})-m_{j}||^{2}$}\\\ 0&\text{otherwise}\end{cases}$

M step. Re-estimate the cluster means using the current assignations
$r_{i,c}$:

$m_{c}\leftarrow\frac{\sum_{i\in\Gamma_{c}}w_{i}\phi(x_{i})}{s_{c}}$

return _$\boldsymbol{r}$_

Algorithm 4 Weighted kernel $k$-means

###  4.2 An spectral solution for the Weighted kernel $k$-means problem

The spectral decomposition of matrix $W^{\frac{1}{2}}KW^{\frac{1}{2}}$ can
provide a different solution for the weighted kernel $k$-means problem. A
standard result in linear algebra [4] provides a global solution to a relaxed
version of problem (24). By allowing $Y$ to be an arbitrary orthonormal
matrix, we can obtain an optimal Y by taking the top $k$ eigenvectors (i.e.,
the eigenvectors associated to the $k$ largest eigenvalues) of the matrix
$W^{\frac{1}{2}}KW^{\frac{1}{2}}$. Each row of the resulting matrix $Y$ is
then interpreted as an "embedding" version of the original data point in a
lower dimensional space with dimension equal to $k$. The typical way to
proceed is to compute a discrete partition of the embedded data points,
usually using the vanilla $k$-means algorithm with the embedded points.

####  4.2.1 The NJW-SC algorithm

The NJW-SC (Ng-Jordan-Weiss Spectral Clustering) algorithm, shown in algorithm
5, considers the weighted kernel $k$-means problem, using the Gaussian kernel

$K_{ij}=\exp{\Big{(}-\frac{\left\lVert
x_{i}-x_{j}\right\rVert^{2}}{2\sigma^{2}}\Big{)}}$ (26)

and setting

$W^{-1}=diag(K1_{n})$ (27)

where $diag(v)$ is a square diagonal matrix with the elements of vector $v$ on
the diagonal, and $1_{n}$ denotes the vector of ones of size $n$. Then, it
proceeds to compute the spectral solution mentioned in the last section: it
computes the top $k$ eigenvectors of the matrix
$W^{\frac{1}{2}}KW^{\frac{1}{2}}$, and use the embedded data points to compute
a discrete partition of the data.

####  4.2.2 NJW-SC complexity

Using a kernel gives the NJW-SC algorithm more flexibility than the vanilla
$k$-means algorithm. However, the cost to pay for this flexibility is the time
complexity of the algorithm. In the first step, the Kernel matrix step, we
need to compute the kernel function for $\frac{n(n+1)}{2}$ pairs of points.
Therefore, the complexity of this step is $O(n^{2}d)$, where $O(d)$ is the
complexity of computing the Gaussian kernel $\kappa_{G}()$. In the second
step, the Weights matrix step, we need to compute the matrix of weights $W$
using equation (27), which implies a matrix-vector multiplication, with
complexity $O(n^{2})$. The third step, the Eigendecomposition step, consist of
two sub-steps. The first one is the matrix multiplication
$W^{\frac{1}{2}}KW^{\frac{1}{2}}$, which can take up to $O(n^{3})$ operations,
but given that $W$ is a diagonal matrix, this sub-step has a complexity of
$O(n^{2})$. The second sub-step, the computation of the top $k$ eigenvectors,
can take up to $O(n^{3})$ operations if computed naively (i.e., computing the
full eigendecomposition). However, as we are interested only in the top $k$
eigenvectors, this sub-step can be computed with more efficient techniques,
such as the power method, that scale like $O(n^{2}k)$ [1]. The complexity of
the last step, the Partition step, depends on the method used to create the
partition. If vanilla $k$-means is used for this purpose, the complexity of
this step is $O(nk^{2}T)$. As we can see, and assuming $k<<d,n$ (which
typically is true in practice), the complexity of the algorithm is
$O(n^{2}d)$, due the computation of the kernel matrix in the first step. This
complexity and other important properties of the NJW-SC algorithm are
summarized in table 3.

Input: Dataset $X$, number of clusters $k$

Output: Learned values for $r_{i,c}$

Kernel matrix step. Compute kernel matrix $K$ for dataset $X$:

$K_{ij}\leftarrow\exp{\Big{(}-\frac{\left\lVert
x_{i}-x_{j}\right\rVert^{2}}{2\sigma^{2}}\Big{)}};$

Weights matrix step. Compute matrix of weights $W$:

$W\leftarrow diag^{-1}(K1_{n});$

Eigendecomposition step. Compute top $k$ eigenvectors:

$Y\leftarrow topEig(W^{\frac{1}{2}}KW^{\frac{1}{2}},k);$

Partition step. Using the embedded points in $Y$, compute a partition
$\boldsymbol{r}$ of the data

return _$\boldsymbol{r}$_

Algorithm 5 NJW-SC Table 3: NJW-SC algorithm | SC
---|---
Number of hyperparameters | 2
Hyperparameters | number of clusters $k$, width of Gaussian kernel $\sigma$
Time complexity | $O(n^{2}d)$
Outlier detection | No
Type of assignment | hard
distance/similarity measure | Gaussian kernel $\kappa_{G}()$ (See Spectral Clustering section for details)

###  4.3 From graphs to Spectral Clustering

In [2], Dhillon, Guan and Kulis present a different formulation of the NJW-SC
algorithm using a graph cut approach. This formulation will be useful for next
section.

####  4.3.1 Normalized graph cuts

Given a graph $G=(\mathcal{V},\mathcal{E},A)$, where $\mathcal{V}$ is a set of
$n$ nodes (or vertices), $\mathcal{E}$ is the set of edges connecting nodes,
and $A\in\mathbb{R}^{n\times n}$ is a non-negative and symmetric edge
similarity matrix, i.e., $A_{ij}$ is a (non-negative) weight expressing the
similarity between node $i$ and node $j$, the $k$-way normalized cut problem
is defined as follows. Suppose $\mathcal{A},\mathcal{B}\subseteq\mathcal{V}$,
and define:

$links(\mathcal{A},\mathcal{B})=\sum_{i\in\mathcal{A},j\in\mathcal{B}}A_{ij}$
(28)

i.e., $links(\mathcal{A},\mathcal{B})$ is the sum of the weights of the edges
that cross from subset $\mathcal{A}$ to subset $\mathcal{B}$, and:

$normlinks(\mathcal{A},\mathcal{B})=\frac{links(\mathcal{A},\mathcal{B})}{links(\mathcal{A},\mathcal{V})}$
(29)

is the normalized link ratio of $\mathcal{A},\mathcal{B}$. The $k$-way
normalized cut problem is to minimize the links that "escape" a cluster
$\mathcal{V}_{c}\subseteq\mathcal{V}$ relative to the total "weight" of the
cluster, where the set of clusters $\mathcal{V}_{c}$ is a partition of
$\mathcal{V}$:

$\displaystyle\min_{\mathcal{V}_{c}}$
$\displaystyle\frac{1}{k}\sum^{k}_{c=1}normlinks(\mathcal{V}_{c},\mathcal{V}\setminus\mathcal{V}_{c})$
(30)

where $\mathcal{V}\setminus\mathcal{V}_{c}$ denotes the set of nodes in
$\mathcal{V}$ that do not belong to cluster $\mathcal{V}_{c}$. We can rewrite
the objective function of problem (30) as follows

$\displaystyle\frac{1}{k}\sum^{k}_{c=1}normlinks(\mathcal{V}_{c},\mathcal{V}\setminus\mathcal{V}_{c})$
$\displaystyle=\frac{1}{k}\sum^{k}_{c=1}\frac{links(\mathcal{V}_{c},\mathcal{V}\setminus\mathcal{V}_{c})}{links(\mathcal{V}_{c},\mathcal{V})}$
$\displaystyle=\frac{1}{k}\sum^{k}_{c=1}\frac{\sum_{i\in\mathcal{V}_{c},j\notin\mathcal{V}_{c}}A_{ij}}{\sum_{i\in\mathcal{V}_{c},j\in\mathcal{V}}A_{ij}}$
$\displaystyle=\frac{1}{k}\sum^{k}_{c=1}\frac{\sum_{i\in\mathcal{V}_{c},j\in\mathcal{V}}A_{ij}-\sum_{i,j\in\mathcal{V}_{c}}A_{ij}}{\sum_{i\in\mathcal{V}_{c},j\in\mathcal{V}}A_{ij}}$
$\displaystyle=\frac{1}{k}\sum^{k}_{c=1}(1-normlinks(\mathcal{V}_{c},\mathcal{V}_{c}))$
$\displaystyle=1-\frac{1}{k}\sum^{k}_{c=1}\frac{\sum_{i,j\in\mathcal{V}_{c}}A_{ij}}{\sum_{i\in\mathcal{V}_{c},j\in\mathcal{V}}A_{ij}}$

We can observe that problem (30) is equivalent to the following optimization
problem

$\displaystyle\max$
$\displaystyle\sum^{k}_{c=1}\frac{\sum_{i,j\in\mathcal{V}_{c}}A_{ij}}{\sum_{i\in\mathcal{V}_{c},j\in\mathcal{V}}A_{ij}}$
(31)

In other words, we are trying to maximize the normalized similarity among
points in the same cluster, where the normalization term is the "weight" of
the cluster, and the similarity between two points is defined by the
similarity matrix $A$. Letting $D\in\mathbb{R}^{n\times n}$ be the diagonal
matrix with node weights $d_{i}$ in the diagonal, where a node weight $d_{i}$
is defined as the sum of weights $A_{ij}$ for edges adjacent to node $i$,
i.e., $d_{i}=\sum^{n}_{j=1}A{ij}$, and $U\in\\{0,1\\}^{n\times k}$ be the
assignment matrix where $U_{ic}=1$ if node $i$ belongs to cluster
$\mathcal{V}_{c}$ and $0$ otherwise, then we can rewrite the objective
function of problem (31) as

$\displaystyle\sum^{k}_{c=1}\frac{\sum_{i,j\in\mathcal{V}_{c}}A_{ij}}{\sum_{i\in\mathcal{V}_{c},j\in\mathcal{V}}A_{ij}}$
$\displaystyle=\sum^{k}_{c=1}\frac{U^{T}_{c}AU_{c}}{U^{T}_{c}DU_{c}}$
$\displaystyle=Tr\Big{(}(U^{T}DU)^{-\frac{1}{2}}U^{T}AU(U^{T}DU)^{-\frac{1}{2}}\Big{)}$

where $U_{c}$ is the $c$-th column of matrix $U$, corresponding to cluster
$\mathcal{V}_{c}$, and, as can be verified,

$(U^{T}DU)^{-\frac{1}{2}}=\begin{bmatrix}\frac{1}{\sqrt{U^{T}_{1}DU_{1}}}&0&...&0\\\
0&\frac{1}{\sqrt{U^{T}_{2}DU_{2}}}&...&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&...&\frac{1}{\sqrt{U^{T}_{k}DU_{k}}}\\\ \end{bmatrix}$ (32)

and

$(U^{T}AU)=\begin{bmatrix}U^{T}_{1}AU_{1}&U^{T}_{1}AU_{2}&...&U^{T}_{1}AU_{k}\\\
U^{T}_{2}AU_{1}&U^{T}_{2}AU_{2}&...&U^{T}_{2}AU_{k}\\\
\vdots&\vdots&\ddots&\vdots\\\
U^{T}_{k}AU_{1}&U^{T}_{k}AU_{2}&...&U^{T}_{k}AU_{k}\\\ \end{bmatrix}$ (33)

Note that, in order to ensure that all points are assigned to exactly one
cluster, any row of matrix $U$ must have exactly one entry with value one, and
the rest equal to zero, i.e., $U1_{k}=1_{n}$. Also, observe that

$(U^{T}U)=\begin{bmatrix}|\mathcal{V}_{1}|&0&...&0\\\
0&|\mathcal{V}_{2}|&...&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&...&|\mathcal{V}_{k}|\\\ \end{bmatrix}$ (34)

We will denote this matrix as $I_{U^{T}U}$. We can then re-cast the
optimization problem (31) as

$\displaystyle\max_{U\in\\{0,1\\}^{n\times k}}$ $\displaystyle
Tr\Big{(}(U^{T}DU)^{-\frac{1}{2}}U^{T}AU(U^{T}DU)^{-\frac{1}{2}}\Big{)}$ (35)
s.t. $\displaystyle U^{T}U=I_{U^{T}U}$ $\displaystyle U1_{k}=1_{n}$

Finally, letting $\tilde{Z}=D^{\frac{1}{2}}U(U^{T}DU)^{-\frac{1}{2}}$ results
in the following optimization problem

$\displaystyle\max_{\tilde{Z}\in\nu_{\tilde{Z}}^{n\times k}}$ $\displaystyle
Tr\Big{(}\tilde{Z}^{T}D^{-\frac{1}{2}}AD^{-\frac{1}{2}}\tilde{Z}\Big{)}$ (36)
s.t. $\displaystyle\tilde{Z}^{T}\tilde{Z}=I$

where
$\nu_{\tilde{Z}}=\\{0\\}\bigcup\left\\{\sqrt{\frac{d_{i}}{U^{T}_{c}DU_{c}}}\right\\},i\in\Gamma_{c}$.
We will denote problem (36) as the matrix version of the $k$-way normalized
cut problem.

Comparing problem (36) with problem (24), we can observe that they are
equivalent problems: matrices $\tilde{Z}$ and $Y$ are analogous, both
orthonormal assignment matrices with elements equal to zero or to the square
root of a normalized weight corresponding to a node ($d_{i}$) or data point
($w_{i}$). Matrices $A$ and $K$ are both similarity matrices: matrix $A$ in a
explicit way and matrix $K$ as a kernel matrix. Finally, and more
specifically, the NJW-SC algorithm sets the matrix $W=diag^{-1}(K1_{n})$, and,
we can rewrite matrix $D$ as

$D=diag(A1_{n})$ (37)

Therefore, we can think of the NJW-SC algorithm as an algorithm that reduces a
data clustering problem to a graph normalized cut problem: it takes the data
points in $X$ and creates a complete graph $G=(\mathcal{V},\mathcal{E},A)$,
where every pair of nodes is connected by a unique edge. This is shown in
figure 9. The graph shown is for illustration purpose only: as mentioned
before, the NJW-SC algorithm creates a fully connected graph.

(a) Banana dataset

(b) Graph created from the banana dataset1

Figure 9: Banana dataset and the corresponding graph created by a SC algorithm
11footnotetext: Image taken from [20]

For every data point in $X$, the NJW-SC algorithm adds a corresponding node in
$\mathcal{V}$, and it sets the similarity matrix $A$ equal to the kernel
matrix $K$ of the data. Finally, it solves the matrix version of the $k$-way
normalized cut problem (36) computing the top $k$ eigenvectors of matrix
$D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ and using them to create a partition of
the nodes in $\mathcal{V}$. This is shown in algorithm 6.

This interesting connection between graph theory and clustering algorithms
will allow us to find a relation between the NJW-SC and the DBSCAN algorithms.

Input: Dataset $X$, number of clusters $k$

Output: $k$ clusters in $X$

Graph step. Create a complete graph $G=(\mathcal{V},\mathcal{E},A)$ such that,
for every point in $X$, there is a corresponding node in $\mathcal{V}$, and
for any pair of nodes in $\mathcal{V}$ there is an edge in $\mathcal{E}$. Set
the similarity matrix $A$ as:

$A_{ij}\leftarrow\exp{\Big{(}-\frac{\left\lVert
x_{i}-x_{j}\right\rVert^{2}}{2\sigma^{2}}\Big{)}};$

Weights matrix step. Compute matrix of weights $D$:

$D\leftarrow diag(A1_{n});$

Eigendecomposition step. Compute top $k$ eigenvectors:

$\tilde{Z}\leftarrow topEig(D^{-\frac{1}{2}}AD^{-\frac{1}{2}},k);$

Partition step. Using the embedded points in $\tilde{Z}$, compute a partition
$\mathcal{V}_{c}$ of the nodes in $\mathcal{V}$

return _$\mathcal{V}_{c}$ , $\forall c\in\\{1,2,...,k\\}$_

Algorithm 6 NJW-SC, Graph version

##  5 DBSCAN

The main difference between DBSCAN and any of the previous algorithms explored
is that DBSCAN is a density based algorithm: it uses two parameters $\epsilon$
and $minPts$ to compute the local density of every data point $x_{i}$ and uses
this information to create a set of cluster from the data. Observe that the
word density is being used loosely: in reality, what DBSCAN computes is the
number of neighbors inside its $\epsilon$-neighborhood
$N_{\epsilon}(x_{i})=\\{x_{j}\in X|\left\lVert
x_{i}-x_{j}\right\rVert\leq\epsilon\\}$. If the number of neighbors in
$N_{\epsilon}(x_{i})$ is equal or larger than $minPts$ (which stands for
minimum number of points), then we label $x_{i}$ as a core point (or dense
point).

###  5.1 SC but faster (and sparser)

Three of the main disadvantages of the NJW-SC algorithm are:

  * •

It needs to find the top $k$ eigenvectors of a $n\times n$ matrix. Despite not
being the term that dominates the complexity of the algorithm, this step adds
$O(n^{2}k)$ computations to the algorithm

  * •

Choice of $k$. Without any prior knowledge, the choice of $k$ can be difficult

  * •

No outlier detection. The algorithm can’t detect outliers. It will assign any
outlier to some cluster, or it will create a size-one cluster with the outlier

We will see how using a different kernel, and introducing a "filtering" step
into the NJW-SC algoritm we can create an algorithm with none of the mentioned
disadvantages.

####  5.1.1 A discontinuous kernel

The Gaussian kernel, used in the NJW-SC algorithm, is a continuous kernel: the
function $\kappa_{G}()$ is continuous w.r.t. its argument $\left\lVert
x_{i}-x_{j}\right\rVert^{2}$. Using this kernel have some advantages: all the
similarity values lies in the $[0,1]$ interval, and the kernel returns $1$ if
and only if $x_{i}=x_{j}$. However, it also has some disadvantages. One of
them is that, no matter how far two points are, the similarity between them is
never $0$. Zero is only achieved in the limit $\left\lVert
x_{i}-x_{j}\right\rVert^{2}\rightarrow\infty$. Using a different kernel we
might overcome this disadvantage. Let $\kappa_{H}()$ be the Heaviside kernel
(based on the Heaviside function), defined as:

$\kappa_{H}(x_{i},x_{j})=\begin{cases}1&\text{if $\epsilon-\left\lVert
x_{i}-x_{j}\right\rVert\geq 0$}\\\ 0&\text{otherwise}\end{cases}$ (38)

From this form, we can see that the Heaviside kernel is discontinuous, as
shown in figure 10. This kernel assigns a similarity value of one to any two
data points if the distance between them is less or equal than $\epsilon$, and
zero otherwise. Using this discontinuous kernel will allows to modify the NJW-
SC algorithm in order to create a faster algorithm.

Figure 10: Gaussian and Heaviside kernels

####  5.1.2 Core points

Modifying algorithm 6 using the Heaviside kernel instead of the Gaussian
kernel is the first modification that we will do to the NJW-SC algorithm. The
second modification will be the following. After the Weights matrix step, we
introduce a filtering step: we will filter the nodes in $G$ using their
weights $d_{i}=D_{ii}$ (in this case, due the use of the Heaviside kernel, the
weight $d_{i}$ is equal to the degree of its corresponding node, plus one).
This filtering process will be simple:

  * •

If $d_{i}=1$ (i.e., if node $i$ is an isolated node), we will label the
corresponding node (and therefore, the corresponding data point $x_{i}$) as an
outlier, removing it from the graph $G$, with all its adjacent edges, and we
will remove the corresponding row and column from the matrix $A$

  * •

If $1<d_{i}\leq minPts$, the corresponding node (and therefore, the
corresponding data point $x_{i}$) will be labeled as unprocessed, removing it
from the graph $G$, with all its adjacent edges, and we will remove the
corresponding row and column from the matrix $A$

  * •

Finally, if $d_{i}>minPts$, we will label the corresponding node (and
therefore, the corresponding data point $x_{i}$) as a core node, preserving it
in the graph $G$, with all its adjacent edges, and preserving the
corresponding row and column in the matrix $A$

We have three different type of data points after the filtering process:
(current) outliers, unprocessed points and core points. We will see that, once
the algorithm has finished, any unprocessed point will be either a part of a
cluster or an outlier. Also, observe that we have reduced the size of matrices
$A$ and $D$ from $n\times n$ to $n_{core}\times n_{core}$, where $n_{core}$ is
the number of core points found in the dataset. From this fact, we can see
that we had improved the complexity of the next step (the Eigendecomposition
step) from $O(n^{2}k)$ to $O(n_{core}^{2}k)$. Also, observe that, although we
haven’t done it explicitly, if two points have a similarity value of zero, the
corresponding edge will be removed (or more precisely, not taken into account)
from graph $G$. Therefore, we can think of (post-filtered) matrix $A$ as the
sum of the adjacency matrix $A_{G}$ of (post-filtered) graph $G$ plus the
identity matrix of appropiate size

$A=A_{G}+I$ (39)

Similarly, we can think of (post-filtered) matrix $D$ as the sum of the degree
matrix $D_{G}$ of (post-filtered) graph $G$ plus the identity matrix of
appropiate size

$D=D_{G}+I$ (40)

There is something important to highlight with respect to this last point.
Since matrix A has been modified, the post-filtered matrix $D$ needs to be
recomputed $D=diag(A1_{n})$. The Laplacian matrix
$\mathcal{L}\in\mathbb{R}^{n_{core}\times n_{core}}$ of (post-filtered) graph
$G$ is defined as:

$\mathcal{L}=D_{G}-A_{G}=D_{G}-A_{G}+(I-I)=D-A$ (41)

where the last equality is due equations (39) and (40). Given the Laplacian
matrix $\mathcal{L}$ and the degree matrix $D_{G}$ of a graph $G$, the
normalized Laplacian of a matrix is defined as:

$\mathcal{L}_{norm}=D_{G}^{-\frac{1}{2}}\mathcal{L}D_{G}^{-\frac{1}{2}}$ (42)

A very well known result in Graph Theory is that eigenvectors related to zero
eigenvalues of the normalized Laplacian indicate connected components in the
graph $G$ [6]. The results still holds if we use the matrix $D$ instead of
$D_{G}$ as can be easily verified. Therefore, given an eigenvector $v_{0}$
related to a zero eigenvalue

$\displaystyle 0$ $\displaystyle=\mathcal{L}_{norm}v_{0}$
$\displaystyle=D_{G}^{-\frac{1}{2}}\mathcal{L}D_{G}^{-\frac{1}{2}}v_{0}$
$\displaystyle=D^{-\frac{1}{2}}\mathcal{L}D^{-\frac{1}{2}}v_{0}$
$\displaystyle=D^{-\frac{1}{2}}(D-A)D^{-\frac{1}{2}}v_{0}$
$\displaystyle=v_{0}-D^{-\frac{1}{2}}AD^{-\frac{1}{2}}v_{0}$

which implies

$D^{-\frac{1}{2}}AD^{-\frac{1}{2}}v_{0}=v_{0}$ (43)

i.e., normalized Laplacian’s eigenvector $v_{0}$ is also an eigenvector of the
matrix $D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$, in this case with eigenvalue equal
to one.

This last observation has an important consequence: the eigenvectors obtained
in the next step (Eigendecomposition step) with associated eigenvalue equal to
one will be indicator vectors for the connected components in (post-filtered)
graph $G$: for any of these eigenvectors, all the non-zero entries corresponds
to points belonging to the same cluster. This property will allows to modify
the Eigendecomposition step in order to remove the dependency to the parameter
$k$.

Just as mentioned in the previous paragraph, in the Eigendecomposition step,
instead of using the $k$ top eigenvectors, we will compute all the
eigenvectors with eigenvalue equal to one. These eigenvectors will be
indicator vectors for the connected components in (post-filtered) $G$, which
in turn will be the clusters formed by the core points. Due this, we can
remove the last step of the algorithm (the Partition step), because we can
read the partition directly from the eigenvectors with eigenvalue equal to
one.

Once we have found the clusters in the set of core points, we need to process
the unprocessed points. Our modified algorithm does this in a simple way: it
will assign each unprocessed point into the same cluster as their closest core
point, if this core point is part of its $\epsilon$-neighborhood. Otherwise,
it will label it as an outlier.

This modified algorithm is shown in algorithm 7, and is equivalent to the
DBSCAN algorithm. However, other versions of the DBSCAN algorithm doesn’t
compute the connected components of $G$ using the eigendecomposition of matrix
$D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$. For instance, in [8], they use instead a
depth-first search (DFS), a more efficient procedure with complexity
$O(n_{core}^{2})$. In algorithm 8, we present a different formulation of the
DBSCAN algorithm, similar to the version presented in [8], and equivalent to
the original formulation of the algorithm in [7]. We prefer to present this
version because it is more similar to algorithm 7.

Input: Dataset $X$, radius $\epsilon$, and minimum number of points $minPts$

Output: Dense-connected clusters and outliers in $X$

Graph step. Create a graph $G=(\mathcal{V},\mathcal{E},A)$ such that, for
every point in $X$, there is a corresponding node in $\mathcal{V}$, and for
any pair of nodes in $\mathcal{V}$ there is an edge in $\mathcal{E}$. Set the
similarity matrix $A$ as:

$A_{ij}\leftarrow\begin{cases}1&\text{if $\epsilon\geq\left\lVert
x_{i}-x_{j}\right\rVert$}\\\ 0&\text{otherwise}\end{cases};$

Weights matrix step. Compute matrix of weights $D$:

$D\leftarrow diag(A1_{n});$

Core points identification step. Identify the core points:

$label_{i}\leftarrow\begin{cases}core&\text{if $d_{i}>minPts$}\\\
unprocessed&\text{if $1<d_{i}\leq minPts$}\\\ outlier&\text{otherwise}\\\
\end{cases}$

If $label_{i}$ is equal to $unprocessed$ or $outlier$, remove node $i$ from
graph $G$, with all its adjacent edges.Make the corresponding changes to
matrices $A$ and $D$

Eigendecomposition step. Compute eigenvectors with eigenvalue equal to one:

$\tilde{Z}\leftarrow EigOne(D^{-\frac{1}{2}}AD^{-\frac{1}{2}});$

The corresponding clusters can be read from $\tilde{Z}$

Outliers identification step. Assign each $unprocessed$ point in $X$ to the
same cluster as their closest core point, if this core point is part of its
$\epsilon$-neighborhood. Otherwise, label it as an outlier:

$label_{i}\leftarrow\begin{cases}core&\text{if closest core point $\in
N_{\epsilon}(x_{i})$}\\\ outlier&\text{otherwise}\end{cases}$

return _Connected components and outliers in $G$_

Algorithm 7 DBSCAN, Spectral version

####  5.1.3 DBSCAN complexity

Identifying the core points in $X$ takes up to $O(n^{2}d)$ operations, because
we need to compute the distance between any pair of points in the dataset.
Constructing the graph in the Graph step takes $O(n_{core}^{2})$ operations,
because that is the maximum possible number of edges. Computing the connected
components of $G$ using DFS instead of the eigendecomposition reduces the
complexity of the Connected components step from $O(n_{core}^{2}k)$ to
$O(n_{core}^{2})$. Finally, and assuming that we have stored in memory the
distances between any pair of points in $X$, the Outliers identification step
takes up to $O(n_{core}(n-n_{core}))$ operations, because we need to find the
closest core point for every unprocessed point. It is evident that the Core
points identification step is the most computational intensive step in
algorithm 8. Therefore, the complexity of DBSCAN is $O(n^{2}d)$. This
complexity and other important properties of DBSCAN are summarized in table 4.

Table 4: DBSCAN algorithm | DBSCAN
---|---
Number of hyperparameters | 2
Hyperparameters | radius $\epsilon$, minimum number of points $minPts$
Time complexity | $O(n^{2}d)$
Outlier detection | Yes
Type of assignment | hard
distance/similarity measure | Heaviside kernel $\kappa_{H}()$ (See DBSCAN section for details)

Input: Dataset $X$, radius $\epsilon$, and minimum number of points $minPts$

Output: Dense-connected clusters and outliers in $X$

Core points identification step. For every $x_{i}$ in $X$, find the points in
its $\epsilon$-neighborhood $N_{\epsilon}(x_{i})$, and identify the core
points with more than $minPts$ neighbors:

$label_{i}\leftarrow\begin{cases}core&\text{if $|N_{\epsilon}(x_{i})|\geq
minPts$}\\\ unprocessed&\text{otherwise}\end{cases}$

Graph step. Create a graph $G=(\mathcal{V},\mathcal{E})$ such that, for every
core point in $X$, there is a corresponding "core"node in $\mathcal{V}$, and
add an edge to $\mathcal{E}$ between this node and all "core" nodes in its
$\epsilon$-neighborhood $N_{\epsilon}(x_{i})$

Connected components step. Find the connected components of graph $G$

Outliers identification step. Assign each non-core point in $X$ to the same
cluster as their closest core point, if this core point is an
$\epsilon$-neighbor. Otherwise, label it as an outlier:

$label_{i}\leftarrow\begin{cases}core&\text{if closest core point $\in
N_{\epsilon}(x_{i})$}\\\ outlier&\text{otherwise}\end{cases}$

return _Connected components and outliers in $G$_

Algorithm 8 DBSCAN, Graph version

###  5.2 Another way to find core points: climbing the hill

One common assumption in Machine Learning is that the data points
$x_{i}\in\mathbb{R}^{d}$ are i.i.d. samples from a unknown density function
$\mathcal{D}$. We can think of this density function as a hypersurface in
$\mathbb{R}^{d+1}$. The peaks of this hypersurface corresponds to the denser
parts of the space $\mathbb{R}^{d}$, i.e., where data points are more
concentrated. This is illustrated in figure 11. In the example shown, data is
drawn from a GMM with three components. Figure 11(a) shows the sampled data
points. Figure 11(b) shows the density contour lines. Figure 11(c) shows the
density function, and a hyperplane representing a level of density, yielding
three clusters.

(a) Data points in two dimensions

(b) Contour plot

(c) Density model with level set

Figure 11: Dataset with its corresponding contour lines and density model

We will use this idea to present a different way to find the core points in a
dataset $X$ sampled from a density function $\mathcal{D}$. The basic idea is
quite simple: every data point $x_{i}$ in $X$ will "climb the hill", trying to
reach it closest peak (i.e., local maxima), as shown in image 12, until the
climbing point reaches a desired level of density (i.e., until it reaches a
desired "altitude"). Because we don’t know how the density function looks
like, we will use the data itself to approximate the shape of the density
function and the desired altitude.

Figure 12: Example of a data point "climbing" the density function

In order to "climb the hill", we will follow a simple procedure proposed by
Yizong in [10]: given a dataset $X$, a data point $x_{i}$ and a radius
$\epsilon$, they first compute the $\epsilon$-neighborhood
$N_{\epsilon}(x_{i})=\\{x_{j}\in X|\left\lVert
x_{i}-x_{j}\right\rVert\leq\epsilon\\}$, and then, using this local
information, they compute the mean of the points in $N_{\epsilon}(x_{i})$.
This mean is the first stop in the path to the peak, and the algorithm
iterates until the point reaches the desired altitude, for every data point. A
detailed explanation on why this simple procedure is guaranteed to climb the
hill is given in [9] and [10]. The main idea is that this procedure estimates
the gradient of the density function $\mathcal{D}$. To see this, consider the
kernel density estimator [1, 9] $f_{\kappa}(x)$:

$f_{\kappa}(x)=\frac{1}{nh^{d}}\sum_{i=1}^{n}\kappa\Bigg{(}\left\lVert\frac{x-x_{i}}{h}\right\rVert^{2}\Bigg{)}$
(44)

where $h^{D}$ is the volume of a hypercube of side $h$ in a $d$ dimensional
space. It can be seen that, as its name suggest, $f_{\kappa}(x)$ is an
estimator of the value of the density $\mathcal{D}$ at the point $x$. This
estimation is done using the sampled data points $x_{i}$. Given the estimator
$f_{\kappa}(x)$, we can compute its gradient w.r.t. $x$

$\displaystyle\nabla f_{\kappa}(x)$
$\displaystyle=\frac{2}{nh^{d+2}}\sum_{i=1}^{n}(x-x_{i})\kappa^{\prime}\Bigg{(}\left\lVert\frac{x-x_{i}}{h}\right\rVert^{2}\Bigg{)}$
$\displaystyle=\frac{2}{nh^{d+2}}\Bigg{[}\sum_{i=1}^{n}\kappa^{\prime}\Bigg{(}\left\lVert\frac{x-x_{i}}{h}\right\rVert^{2}\Bigg{)}\Bigg{]}\Bigg{[}x-\frac{\sum_{i=1}^{n}x_{i}\kappa^{\prime}\big{(}\left\lVert\frac{x-x_{i}}{h}\right\rVert^{2}\big{)}}{\sum_{i=1}^{n}\kappa^{\prime}\big{(}\left\lVert\frac{x-x_{i}}{h}\right\rVert^{2}\big{)}}\Bigg{]}$

where $\kappa^{\prime}()$ is the derivative of $\kappa()$. We can observe that
the first term of $\nabla f_{\kappa}(x)$ is proportional to
$f_{\kappa^{\prime}}(x)$ (where $f_{\kappa^{\prime}}(x)$ is a kernel density
estimator using $\kappa^{\prime}()$ as a kernel instead of $\kappa()$). The
second term is the (negative) mean shift [9]:

$m(x)=x-\frac{\sum_{i=1}^{n}x_{i}\kappa^{\prime}\big{(}\left\lVert\frac{x-x_{i}}{h}\right\rVert^{2}\big{)}}{\sum_{i=1}^{n}\kappa^{\prime}\big{(}\left\lVert\frac{x-x_{i}}{h}\right\rVert^{2}\big{)}}$
(45)

which is the difference between $x$ and a weighted mean of the data points
$x_{i}$. Using the last two equations

$\nabla f_{\kappa}(x)=\frac{2}{h^{2}}f_{\kappa^{\prime}}(x)m(x)$ (46)

yielding

$m(x)=\frac{h^{2}}{2}\frac{\nabla f_{\kappa}(x)}{f_{\kappa^{\prime}}(x)}$ (47)

This last expression shows that the mean shift vector $m(x)$ is proportional
to the gradient of the kernel density estimator $f_{\kappa}(x)$ and therefore,
it is an estimate of the gradient of $\mathcal{D}$ at the point $x$. In order
to climb the hill, we want to move in a direction proportional to $m(x)$:

$\displaystyle x_{i}^{new}$ $\displaystyle=x-m(x)$
$\displaystyle=\frac{\sum_{i=1}^{n}x_{i}\kappa^{\prime}\big{(}\left\lVert\frac{x-x_{i}}{h}\right\rVert^{2}\big{)}}{\sum_{i=1}^{n}\kappa^{\prime}\big{(}\left\lVert\frac{x-x_{i}}{h}\right\rVert^{2}\big{)}}$
$\displaystyle=\frac{\sum_{i=1}^{n}x_{i}\kappa_{H}(x,x_{i})}{\sum_{i=1}^{n}\kappa_{H}(x,x_{i})}$
$\displaystyle=\frac{\sum_{x_{j}\in
N_{\epsilon}(x_{i})}x_{j}}{|N_{\epsilon}(x_{i})|}$

where we have set $h=1$ and used the Heaviside kernel $\kappa_{H}()$ as
$\kappa^{\prime}()$ in the third equality.

Algorithm 9 shows DBSCAN with this modification. Only the first step is
different, the rest of the steps are the same. For each data point $x_{i}$,
the climbing procedure iterates until $|N_{\epsilon}(x_{i})|\geq minPts$,
i.e., until it has reached the desired altitude, or until $x_{i}^{new}=x_{i}$.
Observe that we need to create a copy of the dataset $X$, because we need to
perform the climbing w.r.t. the original (and static) dataset $X$. Also,
observe that the last step (Outliers identification step) is not needed
anymore: after the climbing procedure, every data point in $X_{climb}$ is a
core point because it has climbed until a desired level of density. If some
$x_{i}$ in $X_{climb}$ couldn’t climb to this altitude, it means it is an
outlier, and the Connected components step will find it as a size-one cluster.

Input: Dataset $X$, radius $\epsilon$, and minimum number of points $minPts$

Output: Dense-connected clusters and outliers in $X$

Copy step. Make a copy of the original dataset $X$:

$X_{climb}\leftarrow X$

Core points identification step. For every $x_{i}$ in $X_{climb}$:

while  _$|N_{\epsilon}(x_{i})| <minPts$ or $x_{i}^{new}=x_{i}$_ do

$x^{new}_{i}\leftarrow\frac{\sum_{x_{j}\in
N_{\epsilon}(x_{i})}x_{j}}{|N_{\epsilon}(x_{i})|}$

Graph step. Create a graph $G=(\mathcal{V},\mathcal{E})$ such that, for every
core point in $X_{climb}$, there is a corresponding "core"node in
$\mathcal{V}$, and add an edge to $\mathcal{E}$ between this node and all
nodes in its $\epsilon$-neighborhood $N_{\epsilon}(x_{i})$

Connected components step. Find the connected components of graph $G$

return _Connected components (and outliers) in $G$_

Algorithm 9 DBSCAN, Climbing hill version

We will find how this version of DBSCAN is related to the Mean Shift algorithm
in the next section.

##  6 Mean Shift (MS)

There are different versions of the Mean Shift algorithm. The common point
between them is the climbing procedure, which is very similar to the climbing
procedure mentioned in the previous section. The main difference between them
is the clustering strategy: different versions of the algorithm use different
criteria to create clusters. In this document, the clustering criteria will be
the same criteria as the climbing hill version of DBSCAN, presented in the
previous section.

The main and only difference between the climbing version of DBSCAN presented
in the last section and the Mean Shift algorithm, is that, during the climbing
procedure, points in the Mean Shift algorithm doesn’t stop to climb once a
desired level of density is reached: the climbing points keeps climbing until
they reach their closest local maxima, as shown in algorithm 10. This simple
difference produces (most of the times) very different results to the ones
found by DBSCAN.

Figure 13: In MS, a data point keeps "climbing" until it reachs a "peak"

###  6.1 Mean Shift complexity

Copying the original dataset $X$ in the Copy step takes $O(nd)$ operations.
Identifying the core points in $X$ using the climbing hill procedure can take
up to $O(n^{2}dT)$ operations, because we need to repeatedly (at most $T$
times) compute the distance between $x_{i}$ in $X_{climb}$ and all points in
$X$, for every point $x_{i}$ in $X_{climb}$. Constructing the graph in the
Graph step takes $O(n_{core}^{2})$ operations, because that is the maximum
possible number of edges. Computing the connected components of $G$ in the
Connected components step using DFS can take up to $O(n_{core}^{2})$
operations. Therefore, the complexity of Mean Shift is $O(n^{2}dT)$. This
complexity and other important properties of Mean Shift are summarized in
table 5.

Input: Dataset $X$, radius $\epsilon$

Output: Dense-connected clusters and outliers in $X$

Copy step. Make a copy of the original dataset $X$:

$X_{climb}\leftarrow X$

Core points identification step. For every $x_{i}$ in $X_{climb}$:

while  _$x_{i}$ hasn’t reached a peak or $x_{i}^{new}=x_{i}$_ do

$x^{new}_{i}\leftarrow\frac{\sum_{x_{j}\in
N_{\epsilon}(x_{i})}x_{j}}{|N_{\epsilon}(x_{i})|}$

Graph step. Create a graph $G=(\mathcal{V},\mathcal{E})$ such that, for every
core point in $X_{climb}$, there is a corresponding "core"node in
$\mathcal{V}$, and add an edge to $\mathcal{E}$ between this node and all
nodes in its $\epsilon$-neighborhood $N_{\epsilon}(x_{i})$

Connected components step. Find the connected components of graph $G$

return _Connected components (and outliers) in $G$_

Algorithm 10 Mean Shift Table 5: Mean Shift algorithm | MS
---|---
Number of hyperparameters | 1
Hyperparameters | bandwith (radius) $\epsilon$
Time complexity | $O(n^{2}dT)$
Outlier detection? | Yes
Type of assignment | hard
distance/similarity measure | Heaviside kernel $\kappa_{H}()$ (See Mean Shift section for details)

##  7 Conclusions and future work

This document presents the relationship between different clustering
algorithms. A connection between DBSCAN and Mean Shift is stablished: DBSCAN
can be viewed as a climbing procedure that stops once a data point has reached
a given value of density. Mean shift follows a similar procedure, but its
stopping criteria is different, it stops until a data point has reached its
closest local maxima. Some of the key differences between the algorithms
presented in this document are:

  * •

The distance/similarity measure used. This choice have a huge impact on the
flexibility of the algorithm, i.e, what possible shapes can the algorithm find

  * •

The optimization method. This choice determines the set of possible solutions
(clusterings) that an algorithm can find. For instance, the solutions found by
the EM algorithm and the eigendecomposition method are usually different, and
each optimization methods has its advantages and disadvantages. EM is usually
faster than the eigendecomposition method, but it usually returns a local
optimum. On the other hand, the eigendecomposition method is guaranted to find
the global optimum of the related matrix problem. For example, in section 4,
the eigendecomposition method finds the global optimum of the relaxed matrix
version of the weighted kernel $k$-means problem

These choices have a strong impact in the trade-off between the flexibility
and the time complexity of each algorithm. Table 6 presents a summary of some
of the main features of each algorithm. Figure 14 answers the question in
figure 1:

  * •

Setting all the covariance matrices equal to $\epsilon I$ and considering the
limit $\epsilon\to 0$ take us from GMM to $k$-means

  * •

Adding flexibility to $k$-means via a Gaussian kernel $\kappa_{G}()$ and
introducing a weight $w_{i}$ for every data point yields a Spectral Clustering
algorithm

  * •

Using a different kernel $\kappa_{H}()$ and introducing a filtering step take
us from Spectral Clustering to DBSCAN

  * •

Climbing to the peak instead of stopping at a certain level take us from
DBSCAN to Mean shift

Figure 14: Relationship between GMM, $k$-means, SC, DBSCAN and MS Table 6: Algorithms comparison | GMM |  $k$-means | SC | DBSCAN | MS
---|---|---|---|---|---
Number of hyperparameters | 1 | 1 | 2 | 2 | 1
Hyperparameters | $k$ | $k$ | $k,\sigma$ | $\epsilon,minPts$ | $\epsilon$
Time complexity | $O(nd^{2}kT+d^{2.807}kT)$ | $O(ndkT)$ | $O(n^{2}d)$ | $O(n^{2}d)$ | $O(n^{2}dT)$
Outlier detection | No | No | No | Yes | Yes
Type of assignment | soft | hard | hard | hard | hard
distance/similarity measure | Mahalanobis | Euclidean | Gaussian | Heaviside | Heaviside

###  7.1 Future work

In future work, we would like to investigate what might be the relationship
between some of the algorithms presented and some other clustering algorithms,
like Hierarchical clustering algorithms or some other non-parametric methods,
like Latent Dirichlet Allocation (LDA) [14] or the Chinese Restaurant Process
(CRP) [15]. Similarly, it would be interesting to explore the possible
relationship between some of the algorithms presented and some recently
developed neural-based unsupervised learning algorithms, like the Stacked
Capsule Autoencoders [16], Deep InfoMax (DIM) [17], Invariant Information
Clustering (IIC) [18], or Deep Embedded Clustering (DEC) [19].

## References

  * [1] Christopher M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006.
  * [2] Inderjit S. Dhillon, Yuqiang Guan, and Brian Kulis. Kernel $k$-means, Spectral Clustering and Normalized Cuts. KDD, 2004.
  * [3] Andrew Ng, Michael Jordan and Yair Weiss. On Spectral Clustering: Analysis and an algorithm. NIPS, 2002.
  * [4] Gene H. Golub and Charles F. Van Loan. Matrix Computations. Johns Hopkins University Press, 1989.
  * [5] Bernhard Scholkopf, Christopher J. C. Burges and Alexander J. Smola Advances in Kernel Methods. Support Vector Learning. Springer, 2006.
  * [6] Erich Schubert, Sibylle Hess, and Katharina Morik. The Relationship of DBSCAN to Matrix Factorization and Spectral Clustering. LWDA, 2018.
  * [7] Martin Ester, Hans-Peter Kriegel, Jörg Sander and Xiaowei Xu. A density-based algorithm for discovering clusters in large spatial databases with noise. KDD, 1996.
  * [8] Jennifer Jang and Heinrich Jiang. DBSCAN++: Towards fast and scalable density clustering. ICML, 2019.
  * [9] Peter Meer and Dorin Comaniciu. Mean Shift: A Robust Approach Toward Feature Space Analysis. IEEE Transactions on pattern analysis and machine intelligence, 24(5):603-619, 2002.
  * [10] Yizong Cheng. Mean Shift, Mode Seeking, and Clustering. IEEE transactions on pattern analysis and machine intelligence, 17(8):790-799, 1995.
  * [11] Lawrence Hubert and Phipps Arabie. Comparing partitions. Journal of classification, 2(1):193-218, 1985.
  * [12] Nguyen Xuan Vinh, Julien Epps, and James Bailey. Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance. Journal of Machine Learning Research, 11(Oct):2837-2854, 2010.
  * [13] Roy De Maesschalck, Jouan-Rimbaud Delphine, and Massart Désiré L. The mahalanobis distance. Chemometrics and intelligent laboratory systems, 50(1):1-18, 2000.
  * [14] David M. Blei, Andrew Ng, and Michael Jordan. Latent Dirichlet Allocation. Journal of Machine Learning Research, 3(4–5):993–1022, 2003.
  * [15] J. Pitman. Combinatorial Stochastic Processes. Lecture Notes for St. Flour Summer School, Springer-Verlag, 2002.
  * [16] Adam Kosiorek, et. al. Stacked Capsule Autoencoders. NeurIPS, 2019.
  * [17] R. Devon Hjelm, et. al. Learning deep representations by Mutual Information Estimation and Maximization. ICLR, 2019.
  * [18] Xu Ji, João F. Henriques, and Andrea Vedaldi. Invariant information clustering for unsupervised image classification and segmentation. Proceedings of the IEEE International Conference on Computer Vision, 2019.
  * [19] Junyuan Xie, Ross Girshick, and Ali Farhadi. Unsupervised Deep Embedding for Clustering Analysis. ICML, 2016.
  * [20] Spectral Clustering. (2020, June 6). Retrieved from http://math.ucdenver.edu/~sborgwardt/wiki/index.php/Spectral_clustering