\documentstyle[psbox,ps,11pt,fullpage]{article}
\input epsf

\title{Kernel Clustering in High-Dimensional Feature Spaces}
\author{Matthew Gray}
\begin{document}
\maketitle

\section{Introduction}

 Support Vector Machines (SVMs) are a clever way to do structural risk
 minimization.  The real trick of SVMs is the utilization of a kernel
 function to represent the dot product in a high dimensional feature
 space.  This construction allows an algorithm which can be formulated
 in tears of dot products to be performed in an arbitrarily
 high-dimensional feature space without the computational burden.

\par
Sch\"olkopf, Smola and M\"uller\footnote{Nonlinear Component Analysis
as a Kernel Eigenvalue Problem, TR \#44, Max-Planck-Institut f\"ur
biologische Kybernetic, 1996} show that this technique can be applied
to principal component analysis.  They further suggest that one might
use the same approach for projection pursuit, k-means clustering, k
nearest neighbor, or any other technique that can be cast in terms of
dot products.

\section{Algorithm}

First we define a non-linear mapping $\Phi({\bf x})$ of feature
vectors ${\bf x}$ into our high dimensional feature space, and a
kernel function $k({\bf x}, {\bf y}) = \Phi({\bf x})\cdot\Phi({\bf
y})$. Naturally, the means for our clustering lie in the span of the
$\Phi({\bf x})$ feature vectors.  We will define the coefficients on
these mapped feature vectors as:

\begin{equation}
{\bf m_\nu} = \sum_{j=1}^{M} \gamma_{\nu j}\Phi({\bf x}_j).
\end{equation}


Now, the distance between a mean ${\bf m}_\nu$ and a mapped vector
$\Phi({\bf x})$ can be expressed in terms of dot products:

\begin{equation}
||\Phi({\bf x}) - \sum_{j=1}^M\gamma_{\nu j} \Phi({\bf x}_j)||^2\\
= k({\bf x}, {\bf x}) - 2\sum^{M}_{j=1} \gamma_{\nu j} k({\bf x}, {\bf x}_j) + \sum^{M}_{i,j=1} \gamma_{\nu i}\gamma_{\nu j} k({\bf x}_i, {\bf x}_j)
\label{dist}
\end{equation}

This allows us to determine which cluster is nearest for any given
point.  For training points, the means should be updated according to:

\begin{equation}
{\bf m}_\alpha^{new} = {\bf m}^{old}_\alpha + \zeta(\Phi({\bf x}_{new}) - {\bf m}_\alpha^{old})
\end{equation}

where $\zeta$ is 1 over the number of points closest to ${\bf m}_\alpha$.
\par
Therefore, to update $\gamma$, we do:

\begin{equation}
\gamma^{new}_{\alpha j} = \left\{  
\begin{tabular}{ll}
$\gamma^{old}_{\alpha j}(1 - \zeta)$& for all $j$ other than the current training point\\
&\\
$\zeta + \gamma^{old}_{\alpha j}(1 - \zeta)$& for $j = $ the current training point\\
\end{tabular}
\right.
\end{equation}

After correcting $\gamma$ for each new training point, we recalculate
cluster assignments, and correct $\gamma$ again, until it converges.
\par
To identify which cluster a test point belongs in, we simply calculate
the distance using equation \ref{dist}.
\par

This k-means implementation was tested on a variety of data with
limited success, while it performed well on non-overlapping clusters,
even those that were not linearly seperable, overlapping clusters, as
would be expected, performed poorly.
\par

Modifying the algorithm to use E-M clustering with fixed width
Gaussians is trivial.  The expectation step is calculated as usual,
using Formula \ref{dist} for distance.  In the maximization step, The
updated means in E-M are moved to the the responsibility weighted
centroid of the data.  Therefore, in our representation of the means,
the coefficients $\gamma_{\nu j}$ are equal to the responsibility
cluster $\nu$ has for data point $j$.  We now iterate to the
expectation phase and iterate until the convergence.

\par

The choice of kernel function is the same constraint as it is under
SVMs.  Reasonable choices include $k({\bf x}, {\bf y}) = ({\bf x}
\cdot {\bf y} +1)^n$, which represents a feature space of polynomials,
and $k({\bf x}, {\bf y}) = e^{-\frac{\|{\bf x}-{\bf
y}\|^2}{\sigma^2}}$, which corresponds to an infinite feature space
composed of features of the form $e^{i \omega {\bf x}}$ for all real
values of $\omega$, and generates RBF decision surfaces.

\par

For all of the problems used below the RBF kernel is used, with
manually optimized parameters.

\section{Toy Problem}

An initial toy problem was generated to see the algorithm work.  Two
distributions in 2d were generated and the algorithm was run on the
combined distribution.  Distribution A was a unit variance Gaussian
centered at [0,0].  Distribution B was a mixture of 4 unit variance
Gaussians at [4,4], [5,5], [6,6], and [7,7].
\par

\begin{figure}[htbp]
\begin{center}
\mbox{\epsfysize=3in \epsfbox{/tmp/c1.eps}}
\end{center}
\caption{The first cluster settled on Distribution A}
\label{clust1}
\end{figure}

The algorithm was run with 2 clusters and the densities for each
cluster are shown in Figures \ref{clust1} and \ref{clust2}.  The
overall distribution was zero-centered, so the distributions are
shifted from their original locations in the figures.

\begin{figure}[htbp]
\begin{center}
\mbox{\epsfysize=3in \epsfbox{/tmp/c2.eps}}
\end{center}
\caption{The second cluster settled on Distribution B}
\label{clust2}
\end{figure}

Given this performance with only 2 fixed width Gaussians in the
high-dimensional feature space, the technique seems to have a lot of
potential.

\section{Application}

The algorithm was tested on a set of hand gesture data from the Media
Lab.  The data was composed of two sets (s1 and s1) of three gestures
(point, stop and wave) each.  Each gesture is composed of some number
of frames (usually about 50).  Each frame is used to generate a feature
vector composed of the 2-dimensional coordinates of the hand and the
coefficients of the first 20 eigenimages of the frame.
\par

A preliminary examination of how well the algorithm produces density
estimates in a real problem is promising.  Figure \ref{dens1} shows
the density estimate of a 2-dimensional subset of the features for the
``pointing'' gesture, using a {\em single} high-dimensional Gaussian.

\begin{figure}[htbp]
\begin{center}
\mbox{\epsfysize=4in \epsfbox{/tmp/distrib1.eps}}
\end{center}
\caption{Density estimate on a subset of the features for one frame of one gesture}
\label{dens1}
\end{figure}

Density estimation was performed for three different training/test
data sets.  First, all 20 gestures of a particular type (from sets s1
\& s2) were combined, and a random set of 100 frames were selected.
The clustering was then performed on each gesture type to generate the
density estimates.  All the frames in the dataset were then classified
using a 3-way Bayesian classifier based on these densities.  

\par

Second, 50 samples from the 10 gestures from dataset s1 were taken as
training data and densities were estimated.  These densities were then
used to classify the frames from datasets s1 and s2.

\par

Given these classifications of the individual frames, it is necessary
to combine them somehow.  An initial ad hoc simple majority voting
scheme was tried initially and yielded surprisingly good results.  Due
to this, other combination techniques were not pursued.

\par

The total errors for the first experiment (training on 100 random
samples and testing on the full data set) were
\par
\begin{tabular}{|l|l|l|}
\hline
{\bf Gesture}&{\bf Per frame Accuracy}&{\bf Ad hoc voting accuracy}\\
\hline
\hline
Point&76\%&100\%\\
\hline
Wave&50\%&100\%\\
\hline
Stop&61\%&90\%\\
\hline
\end{tabular}
\par
For training on 50 samples from dataset s1, and testing on the dataset s1:

\begin{tabular}{|l|l|l|}
\hline
{\bf Gesture}&{\bf Per frame Accuracy}&{\bf Ad hoc voting accuracy}\\
\hline
\hline
Point&62\%&100\%\\
\hline
Wave&79\%&100\%\\
\hline
Stop&68\%&90\%\\
\hline
\end{tabular}
\par
and finally, for 50 samples from dataset s1, and testing on the dataset s2:

\begin{tabular}{|l|l|l|}
\hline
{\bf Gesture}&{\bf Per frame Accuracy}&{\bf Ad hoc voting accuracy}\\
\hline
\hline
Point&65\%&100\%\\
\hline
Wave&20\%&0\%\\
\hline
Stop&57\%&100\%\\
\hline
\end{tabular}

\section{Conclusion}

The technique works quite well, but needs some refinement.  Currently,
it is rather slow, and initial estimates for parameters of the kernel function
should be able to be gotten from the training data.
\par

I also want to see whether or not in certain cases it is similar or
equivalent to other existing techniques.  For example, I think the
single Gaussian fit is essentially equivalent to Parzen windows, but
this needs to be determined with rigor.
\par

Further, clusters should be able to be of variable width, and other
applications are worth exploring.
\end{document}
