---
abstract: 'Associating distinct groups of objects (clusters) with contiguous regions of high probability density (high-density clusters), is central to many statistical and machine learning approaches to the classification of unlabelled data. We propose a novel hyperplane classifier for clustering and semi-supervised classification which is motivated by this objective. The proposed *minimum density hyperplane* minimises the integral of the empirical probability density function along it, thereby avoiding intersection with high density clusters. We show that the minimum density and the maximum margin hyperplanes are asymptotically equivalent, thus linking this approach to maximum margin clustering and semi-supervised support vector classifiers. We propose a projection pursuit formulation of the associated optimisation problem which allows us to find minimum density hyperplanes efficiently in practice, and evaluate its performance on a range of benchmark data sets. The proposed approach is found to be very competitive with state of the art methods for clustering and semi-supervised classification.'
author:
- |
    Nicos G. Pavlidis n.pavlidis@lancaster.ac.uk\
    Department of Management Science\
    Lancaster University\
    Lancaster, LA1 4YX, UK David P. Hofmeyr d.hofmeyr@lancaster.ac.uk\
    Department of Mathematics and Statistics\
    Lancaster University\
    Lancaster, LA1 4YF, UK Sotiris K. Tasoulis s.tasoulis@ljmu.ac.uk\
    Department of Applied Mathematics\
    Liverpool John Moores University,\
    Liverpool, L3 3AF, UK
title: Minimum Density Hyperplanes
---

low-density separation, high-density clusters, clustering, semi-supervised classification, projection pursuit

Introduction
============

We study the fundamental learning problem: [*Given a random sample from an unknown probability distribution with no, or partial label information, identify a separating hyperplane that avoids splitting any of the distinct groups (clusters) present in the sample.*]{} We adopt the cluster definition given by @Hartigan1975 [chap. 11], in which a [*high-density cluster*]{} is defined as a maximally connected component of the level set of the probability density function, $p(\x)$, at level $c \geqslant
0$, $$\mathrm{lev}_c p(\x) = \left\{ \x \in \R^d \, \big{\vert}\, p(\x) > c\right\}.$$ An important advantage of this approach over other methods is that it is well founded from a statistical perspective, in the sense that a well-defined population quantity is being estimated.

However, since $p(\x)$ is typically unknown, detecting high-density clusters necessarily involves estimates of this function, and standard approaches to nonparametric density estimation are reliable only in low dimensions. A number of existing *density clustering* algorithms approximate the level sets of the empirical density through a union of spheres around points whose estimated density exceeds a user-defined threshold [@Walther1997; @CuevasFF2000; @CuevasFF2001; @RinaldoW2010]. The choice of this threshold affects both the shape and number of detected clusters, while an appropriate threshold is typically not known in advance. The performance of these methods deteriorates sharply as dimensionality increases, unless the clusters are assumed to be clearly discernible [@RinaldoW2010]. An alternative is to consider the more specific problem of allocating observations to clusters, which shifts the focus to local properties of the density, rather than its global approximation. The central idea underlying such methods is that if a pair of observations belong to the same cluster they must be connected through a path traversing only high-density regions. Graph theory is a natural choice to address this type of problem. [@AzzaliniT2007; @StuetzleN2010] and [@MenardiA2014] have recently proposed algorithms based on this approach. Even these approaches however are limited to problems of low dimensionality by the standards of current applications [@MenardiA2014].

An equivalent formulation of the density clustering problem is to assume that clusters are separated through contiguous regions of low probability density; known as the [*low-density separation*]{} assumption. In both clustering and semi-supervised classification, identifying the hyperplane with the maximum margin is considered a direct implementation of the low-density separation approach. Motivated by the success of support vector machines (SVMs) in classification, maximum margin clustering (MMC) [@XuNLS2004], seeks the maximum margin hyperplane to perform a binary partition (bi-partition) of unlabelled data. MMC can be equivalently viewed as seeking the binary labelling of the data sample that will maximise the margin of an SVM estimated using the assigned labels.

In a plethora of applications data can be collected cheaply and automatically, while labelling observations is a manual task that can be performed for a small proportion of the data only. Semi-supervised classifiers attempt to exploit the abundant unlabelled data to improve the generalisation error over using only the scarce labelled examples. Unlabelled data provide additional information about the marginal density, $p(\x)$, but this is beneficial only insofar as it improves the inference of the class conditional density, $p(\x|y)$. Semi-supervised classification relies on the assumption that a relationship between $p(\x)$ and $p(\x|y)$ exists. The most frequently assumed relationship is that high-density clusters are associated with a single class (cluster assumption), or equivalently that class boundaries pass through low-density regions (low-density separation assumption). The most widely used semi-supervised classifier based on the low-density separation assumption is the semi-supervised support vector machine (S$^3$VM) [@VapnikS1977; @Joachims1999; @ChapelleZ2005]. S$^3$VMs implement the low-density separation assumption by partitioning the data according to the maximum margin hyperplane with respect to both labelled and unlabelled data.

Encouraging theoretical results for semi-supervised classification have been obtained under the cluster assumption. If $p(\x)$ is a mixture of class conditional distributions, @CastelliC1995 [@CastelliC1996] have shown that the generalisation error will be reduced exponentially in the number of labelled examples if the mixture is identifiable. More recently, [@SinghNZ2008] showed that the mixture components can be identified if $p(\x)$ is a mixture of a finite number of smooth density functions, and the separation between mixture components is large. [@Rigollet2007] considers the cluster assumption in a nonparametric setting, that is in terms of density level sets, and shows that the generalisation error of a semi-supervised classifier decreases exponentially given a sufficiently large number of unlabelled data. However, the cluster assumption is difficult to verify with a limited number of labelled examples. Furthermore, the algorithms proposed by @Rigollet2007 and @SinghNZ2008 are difficult to implement efficiently even if the cluster assumption holds. This renders them impractical for real-world problems [@JiYLJH2012].

Although intuitive, the claim that maximising the margin over (labelled and) unlabelled data is equivalent to identifying the hyperplane that goes through regions with the lowest possible probability density has received surprisingly little attention. The work of [@BenDavidLPS2009] is the only attempt we are aware of to theoretically investigate this claim. @BenDavidLPS2009 quantify the notion of a low-density separator by defining the [*density on a hyperplane*]{}, as the integral of the probability density function along the hyperplane. They study the existence of universally consistent algorithms to compute the hyperplane with minimum density. The maximum hard margin classifier is shown to be consistent only in one dimensional problems. In higher dimensions only a soft-margin algorithm is a consistent estimator of the minimum density hyperplane. @BenDavidLPS2009 do not provide an algorithm to compute low density hyperplanes.

This paper introduces a novel approach to clustering and semi-supervised classification which directly identifies low-density hyperplanes in the finite sample setting. In this approach the density on a hyperplane criterion proposed by [@BenDavidLPS2009] is directly minimised with respect to a kernel density estimator that employs isotropic Gaussian kernels. The density on a hyperplane provides a uniform upper bound on the value of the empirical density at points that belong to the hyperplane. This bound is tight and proportional to the density on the hyperplane. Therefore, the smallest upper bound on the value of the empirical density on a hyperplane is achieved by hyperplanes that minimise the density on a hyperplane criterion. An important feature of the proposed approach is that the density on a hyperplane can be evaluated exactly through a one-dimensional kernel density estimator, constructed from the projections of the data sample onto the vector normal to the hyperplane. This renders the computation of minimum density hyperplanes tractable even in high dimensional applications.

We establish a connection between the minimum density hyperplane and the maximum margin hyperplane in the finite sample setting. In particular, as the bandwidth of the kernel density estimator is reduced towards zero, the minimum density hyperplane converges to the maximum margin hyperplane. An intermediate result establishes that there exists a positive bandwidth such that the partition of the data sample induced by the minimum density hyperplane is identical to that of the maximum margin hyperplane.

The remaining paper is organised as follows: The formulation of the minimum density hyperplane problem as well as basic properties are presented in Section \[sec:formulation\]. Section \[sec:MaxMarg\] establishes the connection between minimum density hyperplanes and maximum margin hyperplanes. Section \[sec:Methodology\] discusses the estimation of minimum density hyperplanes and the computational complexity of the resulting algorithm. Experimental results are presented in Section \[sec:Exper\], followed by concluding remarks and future research directions in Section \[sec:concl\].

Problem Formulation {#sec:formulation}
===================

We study the problem of estimating a hyperplane to partition a finite data set, $\X = \{ \x_i\}_{i=1}^n \subset \R^d$, without splitting any of the high-density clusters present. We assume that $\X$ is an i.i.d. sample of a random variable $\mathbf{X}$ on $\R^d$, with unknown probability density function $p:\R^d \to \R^+$. A hyperplane is defined as $H(\v,b) {\mathrel{\mathop:}=}\{ \x \in \R^d \,|\, \v \cdot \x = b\}$, where without loss of generality we restrict attention to hyperplanes with unit normal vector, i.e., those parameterised by $(\v,b) \in \B \times \R$, where $\B = \{ \v \in \R^d \, \big|\, \| \v\|=1\}$. Following @BenDavidLPS2009 we define the [*density on the hyperplane*]{} $H(\v,b)$ as the integral of the probability density function along the hyperplane, $$\label{eq1}
I(\v,b) {\mathrel{\mathop:}=}\int_{ H(\v,b)} p(\x) \dd \x.
$$

We approximate $p(\x)$ through a kernel density estimator with isotropic Gaussian kernels, $$\label{eq3}
\hat{p}(\x \vert \X, h^2 I) = \frac{1}{n (2 \pi h^2)^{d/2} } \sum_{i=1}^n \exp \left\{ -\frac{ \|\x - \x_i \|^2}{2 h^2} \right\}.
$$ This class of kernel density estimators has the useful property that the integral in Equation (\[eq1\]) can be evaluated exactly by projecting $\X$ onto $\v$; constructing a one-dimensional density estimator with Gaussian kernels and bandwidth $h$; and evaluating the density at $b$, $$\begin{aligned}
\label{eq4}
\hat I(\v,b| \X, h^2 I) & {\mathrel{\mathop:}=}\int_{H(\v, b)} \hat{p} \left(\x | \X, h^2 I \right) \dd \x, \nonumber \\
& = \frac{1}{n\sqrt{2\pi h^2 }} \sum_{i=1}^n  \exp \left\{- \frac{ (b - \v \cdot \x_i)^2 }{2 h^2} \right\}
=  \hat{p} \left(b \, | \, \{ \v\cdot \x_i \}_{i=1}^n, h^2 \right).
$$ The univariate kernel estimator $\hat{p} \left(\cdot \, | \, \{ \v\cdot \x_i \}_{i=1}^n, h^2 \right)$ approximates the [*projected density on $\v$*]{}, that is, the density function of the random variable, $X_\v = \mathbf{X} \cdot \v$. Henceforth we use $\hat{I}(\v,b)$ to approximate $I(\v,b)$. To simplify terminology we refer to $\hat I(\v,b)$ as the [*density on $H(\v,b)$*]{}, or the [*density integral on $H(\v,b)$*]{}, rather than the empirical density, or the empirical density integral, respectively. For notational convenience we write $\hat{I}(\v, b)$ for $\hat{I}(\v,b| \X, h^2 I)$, where $\X$ and $h$ are apparent from context.

The following Lemma, adapted from [@TasoulisTP2010 Lemma 3], shows that $\hat{I}(\v,b)$ provides an upper bound for the maximum value of the empirical density at any point that belongs to the hyperplane.

\[lem:levsetbound\] Let $\X = \{\x_i\}_{i=1}^n \subset \R^d$, and $\hat{p}(\x|\X,h^2 I)$ be a kernel density estimator with isotropic Gaussian kernels. Then, for any $(\v,b) \in \B\times \R$, $$\label{lm:Tasoulis}
\max_{\x \in H(\v,b)} \hat{p} (\x| \X,h^2 I) \leqslant \left(2 \pi h^2 \right)^{\frac{1-d}{2}} \hat{I}(\v, b),
\quad \textrm{for all } \x \in H(\v, b).
$$

This lemma shows that a hyperplane, $H(\v,b)$, cannot intersect level sets of the empirical density with level higher than $\left(2 \pi h^2
\right)^{\frac{1-d}{2}} \hat{I}(\v, b)$. The proof of the lemma relies on the fact that projection contracts distances, and follows from simple algebra. In Equation (\[lm:Tasoulis\]) equality holds if and only if there exists $\x \in H(\v,b)$ and $\mathbf{c} \in \R^n$ such that all $\x_i \in
\X$, can be written as $\x_i = \x + c_i \v$. It is therefore not possible to obtain a uniform upper bound on the value of the empirical density at points that belong to $H(\v,b)$ that is lower than $\left(2 \pi h^2 \right)^{\frac{1-d}{2}} \hat{I}(\v, b)$ using only one-dimensional projections. Since the upper bound of Lemma \[lem:levsetbound\] is tight and proportional to $\hat{I}(\v, b)$, minimising the density on the hyperplane leads to the lowest upper bound on the maximum value of the empirical density along the hyperplane separator.

To obtain hyperplane separators that are meaningful for clustering and semi-supervised classification, it is necessary to constrain the set of feasible solutions, because the density on a hyperplane can be made arbitrarily low by considering a hyperplane that intersects only the tail of the density. In other words, for any $\v$, $\hat{I}(\v,b)$ can be made arbitrarily low for sufficiently large $|b|$. In both problems the constraints restrict the feasible set to a subset of the hyperplanes that intersect the interior of the convex hull of $\X$. In detail, let $\conv\X$ denote the convex hull of $\X$, and assume $\mbox{Int}(\conv\X)\neq \emptyset$. Define $C$ to be the set of hyperplanes that intersect $\mbox{Int}(\conv\X)$, $$\label{eq:domain}
C =\left\{ H(\v, b) \, \Big{\vert} \, (\v, b) \in \B \times \R, \, \exists \mathbf{z} \in \mbox{Int}(\conv
\X) \mbox{ s.t. } \v \cdot \mathbf{z} = b\right\}.
$$ Then denote by $F$ the set of feasible hyperplanes, where $F \subset C$. We define the [*minimum density hyperplane*]{} (MDH), $H(\vs,\bs) \in F$ to satisfy, $$\label{eq2}
\hat{I}(\vs,\bs) = \min_{(\v, b) \vert H(\v, b) \in F} \hat{I}(\v,b).
$$ In the following subsections we discuss the specific formulations for clustering and semi-supervised classification in turn.

Clustering {#ssec:UL}
----------

Since high-density clusters are formed around the modes of $p(\x)$, the convex hull of these modes would be a natural choice to define the set of feasible hyperplanes. Unfortunately, this convex hull is unknown and difficult to estimate. We instead propose to constrain the distance of hyperplanes to the origin, $b$. Such a constraint is inevitable as for any $\v \in \B$, $\hat{I}(\v,b)$ can become arbitrarily close to zero for sufficiently large $|b|$. Obviously, such hyperplanes are inappropriate for the purposes of bi-partitioning as they assign all the data to the same partition. Rather than fixing $b$ to a constant, we constrain it in the interval, $$\label{constr1}
F(\v) = \left[ \mu_{\v} - \alpha \sigma_\v, \mu_\v + \alpha \sigma_\v \right],
$$ where $\mu_\v$ and $\sigma_\v$ denote the mean and standard deviation, respectively, of the projections $\{\v \cdot \x_i\}_{i=1}^n$. The parameter $\alpha \geqslant 0$, controls the width of the interval, and has a probabilistic interpretation from Chebyshev’s inequality. Smaller values of $\alpha$ favour more balanced partitions of the data at the risk of excluding low density hyperplanes that separate clusters more effectively. On the other hand, increasing $\alpha$ increases the risk of separating out only a few outlying observations. We discuss in detail how to set this parameter in the experimental results section. If $\mathrm{Int}( \conv \X) \neq \emptyset$, then there exists $\alpha>0$ such that the set of feasible hyperplanes for clustering, $F_{\mathrm{CL}}$, satisfies, $$\label{eq:feasible}
F_{\mathrm{CL}} = \left\{ H(\v,b) \, \Big{\vert} \, (\v,b) \in \B \times \R, \, b \in F(\v) \right\} \subset C,
$$ where $C$ is the set of hyperplanes that intersect $\mathrm{Int}(\conv \X)$, as defined in Equation (\[eq:domain\]).

The minimum density hyperplane for clustering is the solution to the following constrained optimisation problem,

\[eq:ULconstr1\] $$\begin{aligned}
& \min_{(\v, b) \in  \B \times \R} \  \hat I(\v, b),  \\
\textrm{subject to:} \ \ &  b - \mu_{\v} + \alpha \sigma_\v \geqslant 0, \\
& \mu_{\v} + \alpha \sigma_\v - b \geqslant 0.
$$

Since the objective function and the constraints are continuously differentiable, MDHs can be estimated through constrained optimisation methods like sequential quadratic programming (SQP). Unfortunately the problem of local minima due to the nonconvexity of the objective function seriously hinders the effectiveness of this approach.

To mitigate this we propose a parameterised optimisation formulation, which gives rise to a projection pursuit approach. Projection pursuit methods optimise a measure of “interestingness” of a linear projection of a data sample, known as the projection index. For our problem the natural choice of projection index for $\v$ is the minimum value of the projected density within the feasible region, $\min_{b \in F(\v)} \hat{I}(\v,b)$. This index gives the minimum density integral of feasible hyperplanes with normal vector $\v$. To ensure the differentiability of the projection index we incorporate a penalty term into the objective function. We define the penalised density integral as, $$\begin{aligned}
\fUL(\v,b) & = \hat{I}(\v,b) + \frac{L}{\eta^\epsilon} \max \left\{ 0, \mu_\v  -\alpha \sigma_\v - b, b - \mu_\v  -\alpha \sigma_\v \right\}^{1+\epsilon},\label{eq:PenalisedUL} 
$$ where, $L = \left(e^{1/2}h^2 \sqrt{2 \pi}\right)^{-1} \geqslant \sup_{b \in \R}
\left| \frac{\partial \hat{I}(\v,b)}{\partial{b}} \right|$, $\epsilon \in (0,1)$ is a constant term that ensures that the penalty function is everywhere continuously differentiable, and $\eta \in
(0,1)$. Other penalty functions are possible, but we only consider the above due to its simplicity, and the fact that its parameters offer a direct interpretation: $L$ in terms of the derivative of the projected density on $\v$; and $\eta$ in terms of the desired accuracy of the minimisers of $\fUL(\v,b)$ relative to the minimisers of Equation (\[eq:ULconstr1\]), as discussed in the following proposition.

\[prop:minim\]

For $\v \in \B$, define, the set of minimisers, $$\begin{aligned}
    B(\v) & =  \operatorname*{arg\,min}_{b \in F(\v)} \hat{I}(\v,b), \\
B_C(\v) & = \operatorname*{arg\,min}_{b \in \R} \fUL(\v,b)
$$ For every $\bs \in B(\v)$ there exists $\bs_C \in B_C(\v)$ such that $| \bs  - \bs_C| \leqslant \eta$. Moreover, there are no minimisers of $\fUL(\v,b)$ outside the interval $[\mu_v - \alpha \sigma_\v -\eta,  \mu_v + \alpha \sigma_\v + \eta]$, $$B_C(\v) \cap \R \backslash [\mu_v - \alpha \sigma_\v -\eta, \mu_v + \alpha \sigma_\v +\eta] = \emptyset.
$$

Any minimiser in the interior of the feasible region, $\bs \in  B(\v) \cap \mbox{Int}(F(\v))$, also minimises the penalised function, since $\fUL(\v,b)= \hat{I} (\v, b)$ for all $b \in \mbox{Int}(F(\v))$, hence $\bs \in B_C(\v)$.

Next we consider the case when either or both of the boundary points of $F(\v)$, $b^- = \mu_v - \alpha \sigma_\v$ and $b^+ = \mu_v + \alpha \sigma_\v$, are contained in $B(\v)$. It suffices to show that, $\fUL(\v,b) > \hat{I}(\v, b^-)$ for all $b < b^- - \eta$, and $\fUL(\v,b) > \hat{I}(\v, b^+)$ for all $b > b^+ + \eta$. We discuss only the case $b > b^+ + \eta$ as the treatment of $b<b^- - \eta$ is identical. Assume that $\hat{I}(\v, b) < \hat{I}(\v, b^+)$ (since in the opposite case the result follows immediately: $\fUL(\v,b) > \hat{I}(\v, b) > \hat{I}(\v, b^+)$). From the mean value theorem there exists $\xi \in (b^+,b)$ such that, $$\begin{aligned}
\hat{I}(\v, b^+) & = \hat{I}(\v, b) - (b - b^+) \left.\frac{ \partial{ \hat{I}(\v,b)} }{\partial{b}} \right|_{b=\xi} \\
& \leqslant \hat{I}(\v,b) + (b - b^+) L \\
& < \hat{I}(\v, b) + \frac{L (b - b^+)^{1+\epsilon}}{\eta^\epsilon} = \fUL(\v,b).
$$ In the above we used the following facts: $\left.\frac{ \partial{ \hat{I}(\v,b)} }{\partial{b}} \right|_{b=\xi} <0$, $L \geqslant \sup_{b \in \R} \left| \frac{ \partial{ \hat{I}(\v,b)} }{\partial{b}} \right|$, and $\frac{b - b^+}{\eta}>1$.

We define the projection index for the clustering problem as the minimum of the penalised density integral, $$\label{eq:PrIndexUL}
\phiUL(\v) = \min_{b \in \R} \fUL(\v,b).
$$ Since the optimisation problem of Equation (\[eq:PrIndexUL\]) is one-dimensional it is simple to compute the set of global minimisers $B_C(\v)$. As we discuss in Section \[sec:Methodology\], this is necessary to compute directional derivatives of the projection index, as well as, to determine whether $\phiUL$ is differentiable. We call the optimisation of $\phiUL$, [*minimum density projection pursuit*]{} (MDP$^2$). For each $\v$, considers only the optimal choice of $b$. This enables it to avoid local minima of $\hat{I}(\v, \cdot)$. Most importantly is able to accommodate a discontinuous change in the location of the global minimiser(s), $\operatorname*{arg\,min}_{b \in\R} \fUL(\v,b)$, as $\v$ changes. Neither of the above can be achieved when the optimisation is jointly over $(\v,b)$ as in the original constrained optimisation problem, Equation (\[eq:ULconstr1\]). The projection index $\phiUL$ is continuous, but it is not guaranteed to be everywhere differentiable when $B_C(\v)$ is not a singleton. The resulting optimisation problem is therefore nonsmooth and nonconvex.

To illustrate the effectiveness of to estimate MDHs, we compare this approach with a direct optimisation of the constrained problem given in Equation (\[eq:ULconstr1\]) using SQP. To enable visualisation we consider the two-dimensional S1 data set [@FrantiV2006], constructed by sampling from a Gaussian mixture distribution with fifteen components, where each component corresponds to a cluster.

Figure \[fig:SQPvsMDP\] depicts the MDHs obtained over 100 random initialisations of SQP and MDP$^2$. It is evident that SQP frequently yields hyperplanes that intersect regions with high probability density thus splitting clusters. As SQP always converged in these experiments the poor performance is solely due to convergence to local minima. In contrast, converges to three different solutions over the 100 experiments, all of which induce high quality partitions, and none intersects a high-density cluster.

![Projection index for S1 data set and solutions obtained through SQP and MDP$^2$[]{data-label="fig:hist"}](histogram2_SQPvsGS.pdf)

In polar coordinates any $\v \in \mathcal{S}^1$ can be parameterised through a single projection angle. Using this parameterisation, the upper plot of Figure \[fig:hist\] depicts the value of the projection index, $\phiUL(\v(\theta))$, for $\theta \in [0,\pi]$. The lower plot of the figure provides histograms of the distribution of the solutions (locally optimal projection angles) obtained over the 100 experiments with SQP (grey) and MDP$^2$ (white). The figure shows that $\phiUL(\v)$ is continuous but not everywhere differentiable. The solution most frequently obtained through MDP$^2$ corresponds to the global optimum, while the only other two solutions identified are the local minimisers with the next two lowest function values. In contrast SQP converges to a much wider range of solutions. Note that this method is not guaranteed to identify the optimal value of $b$ for any $\v(\theta)$ and this indeed occurs in this example. Therefore the value of $\phiUL(\v)$ is a lower bound for the function values of the minimisers identified through SQP.

Semi-Supervised Classification {#ssec:SSL}
------------------------------

In semi-supervised classification labels are available for a subset of the data sample. The resulting classifier needs to predict as accurately as possible the labelled examples, while avoiding intersection with high-density regions of the empirical density. The MDH formulation can readily accommodate partially labelled data by incorporating the linear constraints associated with the labelled data into the clustering formulation. Without loss of generality assume that the first $\ell$ examples are labelled by $\y = (y_1,\ldots,y_\ell)^{\top} \in \{-1,1\}^\ell$. The MDH for semi-supervised classification is the solution to the problem,

\[eq:SSconstr1\] $$\begin{aligned}
& \min_{(\v, b) \in  \B \times \R} \  \hat I(\v, b), \label{eq:SSLobj1} \\
\textrm{subject to:} & \quad  \ y_i(\v \cdot \x_i - b) \geqslant 0, \
    \forall i=1,\ldots,\ell, \label{eq:SSLconstr} \\
    &\quad b - \mu_{\v} + \alpha \sigma_\v \geqslant 0,
    \label{eq:SSLConstr2}\\
&\quad \mu_\v + \alpha \sigma_\v - b \geqslant 0,\label{eq:SSLConstr3}
$$

where $\hat I(\v, b), \mu_\v$, and $\sigma_\v$ are computed over the entire data set. If the labelled examples are linearly separable the constraints in Equation (\[eq:SSconstr1\]) define a nonempty feasible set of hyperplanes, $$F_{\mathrm{LB}} = \left\{ H(\v,b) \,|\, (\v,b) \in \B\times\R, b \in F(\v), y_i(\v \cdot \x_i - b) \geqslant 0, \  \forall i \in \{1,\ldots,l\} \right\}\subset C.  
$$ Equations (\[eq:SSLConstr2\]) and (\[eq:SSLConstr3\]) act as a [*balancing constraint*]{} which discourages MDHs that classify the vast majority of unlabelled data to a single class. Balancing constraints are included in the estimation of §3VMs for the same reason [@Joachims1999; @ChapelleZ2005].

As in the case of clustering, the direct minimisation of Equation (\[eq:SSconstr1\]) frequently leads to locally optimal solutions. To mitigate this we again propose a projection pursuit formulation. We define the penalised density integral for semi-supervised classification as, $$\begin{aligned}
\label{eq:PenalisedSSL}
\fSL(\v,b) = \fUL(\v,b) + \gamma \sum_{i=1}^l \max \left\{0, -y_i (\v \cdot \x_i - b)\right\}^{1+\epsilon}
$$ where, $\gamma>0$ is a user-defined constant, which controls the trade-off between reducing the density on the hyperplane, and misclassifying the labelled examples. The projection index is then defined as the minimum of the penalised density integral, $$\begin{aligned}
\phiSSL(\v) & = \min_{b \in \R} \fSL(\v,b) .\label{eq:PrIndexSSL} 
$$

Connection to Maximum Margin Hyperplanes {#sec:MaxMarg}
========================================

In this section we discuss the connection between MDHs and maximum (hard) margin hyperplane separators. The margin of a hyperplane $H(\mathbf{v}, b)$ with respect to a data set $\X$ is defined as the minimum Euclidean distance between the hyperplane and its nearest datum, $$\mbox{margin}\, H(\mathbf{v}, b) = \min_{\x \in \X}\vert \mathbf{v}\cdot \x - b \vert.
$$ The points whose distance to the hyperplane $H(\v,b)$ is equal to the margin of the hyperplane, that is, $
\operatorname*{arg\,min}_{\x \in \X} \vert \mathbf{v}\cdot \x - b \vert,
$ are called the [*support points*]{} of $H(\v,b)$. Let $F$ denote the set of feasible hyperplanes; then the [*maximum margin hyperplane*]{} (MMH), $H(\v^m,b^m) \in F$ satisfies, $$\margin\, H(\v^m,b^m) = \max_{(\v, b) \vert H(\v, b)\in F} \margin\, H(\v,b).
$$

The main result of this section is Theorem \[thm:convergence\], which states that as the bandwidth parameter, $h$, is reduced to zero the MDH converges to the MMH. An intermediate result, Lemma \[lm:samedivision\], shows that there exists a positive bandwidth, $h'>0$ such that, for all $h \in (0,h')$, the partition of the data set induced by the MDH is identical to that of the MMH.

We first discuss some assumptions which allow us to present the theoretical results of this section. As before we assume a fixed and finite data set $\X \subset \R^d$, and approximate its (assumed) underlying probability density function via a kernel density estimator using Gaussian kernels with isotropic bandwidth matrix $h^2 I$. We assume that the interior of the convex hull of the data, $\mbox{Int}(\conv \X)$, is non-empty, and define $C$ as the set of hyperplanes that intersect $\mbox{Int}(\conv \ \X)$, as in Equation (\[eq:domain\]). The set of feasible hyperplanes, $F$, for either clustering or the semi-supervised classification satisfies $F \subset C$. By construction every $H(\mathbf{v}, b) \in F$ defines a hyperplane which partitions $\X$ into two non-empty subsets. Observe that if for each $\v \in \B$ the set $\{b \in \R \,\vert\, H(\v, b) \in
F\}$ is compact, then by the compactness of $\B$ a maximum margin hyperplane in $F$ exists. For both the clustering and semi-supervised classification problems this compactness holds by construction.

For any $h>0$, let $(\v^\star_h, b^\star_h) \in \B\times\R$ parameterise a hyperplane which achieves the minimal density integral over all hyperplanes in $F$, for bandwidth matrix $h^2I$. That is, $$\hat{I}(\v^\star_h, b^\star_h) = \min_{(\v, b) \vert H(\v, b)\in F} \hat I(\v, b).
$$ Following the approach of @TongK2000 we first show that as the bandwidth, $h$, is reduced towards zero, the density on a hyperplane is dominated by its nearest point. This is achieved by establishing that for all sufficiently small values of $h$, a hyperplane with non-zero margin has lower density integral than any other hyperplane with smaller margin.

\[lm:bddmargin\] Take $H(\v,b) \in F$ with non-zero margin and $0< \delta < \margin\, H(\v,b): = M_{\v, b}$. Then $\exists h' >0$ such that $h \in (0,h')$ and $M_{\mathbf{w}, c}:=\margin\, H(\mathbf{w}, c) \leqslant M_{\v, b}-\delta$ implies $\hat{I}(\v,b) < \hat{I}(\w,c)$.

Using Equation (\[eq4\]) it is easy to see that, $$\begin{aligned}
\hat{I}(\v,b) & \leqslant \frac{1}{ h \sqrt{2 \pi} } \exp \left\{ - \frac{M_{\v, b}^2}{2 h^2} \right\}, \\
\inf \left\{ \hat I(\w,c) \, | \, M_{\w, c} \leqslant M_{\v, b} - \delta \right\} & \geqslant
    \frac{1}{n  h \sqrt{2 \pi} } \exp  \left\{ - \frac{(M_{\v, b}-\delta)^2}{2 h^2} \right\}.
$$ Therefore, $$\begin{aligned}
0 \leqslant \lim_{h \to 0^+} \frac{ \hat I(\v,b) }{\inf \left\{ \hat I(\w,c) \, | \, M_{\w, c} \leqslant M_{\v, b} - \delta \right\}}
& \leqslant \lim_{h \to 0^+} \frac{ n \exp \left\{ - \frac{M_{\v, b}^2}{2 h^2} \right\} }{ \exp \left\{ -\frac{(M_{\v, b}-\delta)^2}{2 h^2} \right\} } = 0.
$$

Therefore, $\exists h'>0$ such that $h \in (0,h^\prime) \Rightarrow \frac{\hat I(\v, b)}{\inf\left\{\hat I(\w, c) \big \vert M_{\w, c} \leqslant M_{\v, b}-\delta \right\}}<1$.

An immediate corollary of Lemma \[lm:bddmargin\] is that as $h$ tends to zero the margin of the MDH tends to the maximum margin. However, this does not necessarily ensure the stronger result that the sequence of MDHs converges to the MMH. To establish this we require two technical results, which describe some algebraic properties of the MMH, and are provided as part of the proof of Theorem \[thm:convergence\] which is given in Appendix \[app:thm:convergence\].

The next lemma uses the previous result to show that there exists a positive bandwidth, $h^\prime>0$, such that an MDH estimated using $h\in (0,h^\prime)$ induces the same partition of $\X$ as the MMH. The result assumes that the MMH is unique. Notice that if $\X$ is a sample of realisations of a continuous random variable then this uniqueness holds with probability 1.

\[lm:samedivision\]

Suppose there is a unique hyperplane in $F$ with maximum margin, which can be parameterised by $(\v^m, b^m) \in \B \times \R$. Then $\exists h'>0$ s.t. $h \in (0,h^\prime) \Rightarrow H(\v^\star_h, b^\star_h)$ induces the same partition of $\X$ as $H(\v^m, b^m)$.

Let $M = \margin\, H(\v^m, b^m)$, and let $P$ be the collection of hyperplanes that induce the same partition of $\X$ as that induced by $H(\v^m, b^m)$. Since $\X$ is finite and $H(\v^m,b^m)$ is unique, $\exists \delta > 0$ s.t. $H(w, c) \notin P
\Rightarrow \margin\, H(w, c) \leqslant M-\delta$. By Lemma \[lm:bddmargin\], $\exists h'>0$ s.t., $$h\in (0,h^\prime)
\Rightarrow H(\v^{\star}_h, b^{\star}_h) \notin \left\{H(\w, c) \, \vert \, \margin\, H(\w, c)
\leqslant M-\delta\right\},$$ therefore $H(\v^{\star}_h, b^{\star}_h) \in P$.

The next theorem is the main result of this section, and states that the MDH converges to the MMH as the bandwidth parameter is reduced to zero. Notice that by the non-unique representation of hyperplanes, the maximum margin hyperplane has two parameterisations in $C$, namely $(\mathbf{v}^m,
b^m)$ and $(-\mathbf{v}^m, -b^m)$. Convergence to the maximum margin hyperplane is therefore equivalent to showing that, $$\min\{\|(\mathbf{v}_h^\star, b_h^\star)-(\mathbf{v}^m, b^m)\|, \|(\mathbf{v}_h^\star, b_h^\star)+(\mathbf{v}^m, b^m)\|\} \to 0 \mbox{ as } h \to 0^+.
$$

\[thm:convergence\]

Suppose there is a unique hyperplane in $F$ with maximum margin, which can be parameterised by $(\v^m, b^m) \in \B \times \R$. Then, $$\lim_{h \to 0^+} \min \left\{\|(\v^\star_h, b^\star_h) - (\v^m, b^m)\|,
    \|(\v^\star_h, b^\star_h) + (\v^m, b^m)\| \right\} = 0.$$

The set $F$ used in Theorem \[thm:convergence\] is generic so it can capture the constraints associated with both clustering and semi-supervised classification, Equations (\[eq:ULconstr1\]) and (\[eq:SSconstr1\]) respectively. In the case of semi-supervised classification we must also assume that the labelled data are linearly separable. Theorem \[thm:convergence\] is not directly applicable to the formulations as in this case the function being minimised is not the density on a hyperplane. The next two subsections establish this result for the formulation of the clustering and semi-supervised classification problem.

for Clustering
--------------

We have shown that for the constrained optimisation formulation the MDH converges to the MMH within the feasible set, $F_{\mathrm{CL}} \subset C$. In addition, for a fixed $\v$, Proposition \[prop:minim\] bounds the distance between minimisers of the penalised function $\fUL$, $\operatorname*{arg\,min}_{b \in \R} \fUL(\v,b)$, and the optimal $b$ of the constrained problem, $\operatorname*{arg\,min}_{b \in F(\v)} \hat{I}(\v,b)$. Combining these we can show that the optimal solution to the penalised MDP$^2$ formulation converges to the maximum margin hyperplane in $F_{\mathrm{CL}}$, provided the parameters within the penalty term suitably depend on the bandwidth parameter, $h$. While the general case can be shown, for ease of exposition we make the simplifying assumption that the maximum margin hyperplane is strictly feasible, i.e., if $(\v^m, b^m)$ parameterises the maximum margin hyperplane then $b^m \in (\mu_{\v^m}-\alpha \sigma_{\v^m},
\mu_{\v^m}+\alpha \sigma_{\v^m})$.

For $h, \eta, L > 0$ define $(\v^\star_{h, \eta, L}, b^\star_{h, \eta, L})$ to be any global minimiser of $f_{\mathrm{CL}}$, i.e., $$f_{\mathrm{CL}}(\v^\star_{h, \eta, L}, b^\star_{h, \eta, L}) = \min_{(\v, b) \in\B\times\R} f_{\mathrm{CL}}(\v, b).$$

\[lem:ULConv\] Suppose there is a unique hyperplane in $F_{\mathrm{CL}}$ with maximum margin, which can be parameterised by $(\v^m, b^m) \in \B \times \R$. Suppose further that $b^m \in (\mu_{\v^m}-\alpha \sigma_{\v^m}, \mu_{\v^m}+\alpha \sigma_{\v^m})$. For $h>0$, let $L(h) = (e^{1/2}h^2\sqrt{2\pi})^{-1}$, and $0 < \eta(h) \leqslant h$. Then, $$\lim_{h \to 0^+} \min\{\|(\v^\star_{h, \eta(h), L(h)}, b^\star_{h, \eta(h), L(h)}) - (\v^m,b^m)\|, \|(\v^\star_{h, \eta(h), L(h)}, b^\star_{h, \eta(h), L(h)}) + (\v^m, b^m)\|\} = 0.$$

Let $M = \mbox{margin}H(\v^m, b^m)$ and as in the proof of Lemma \[lm:samedivision\], let $\delta > 0$ be such that any hyperplane inducing a different partition from $H(\v^m, b^m)$ has margin at most $M-\delta$. Consider the set $F^\delta_{\mathrm{CL}}: = \{(\v, b) \in \B \times \R \vert b \in \mathbb{B}_{\delta/2}(F(\v))\}$, where we used the notation $\mathbb{B}_{\delta/2}(F(\v))$ to denote the neighbourhood of $F(\v)$ given by $\{r \in \R \vert d(r,F(\v))<\delta/2\}$. The set $F^\delta_{\mathrm{CL}}$ increases the feasible set of hyperplanes by allowing $b$ to range in $b \in \mathbb{B}_{\delta/2}(F(\v))$. For any fixed $\v$, the maximum margin of all hyperplanes with normal vector $\v$ can increase by at most $\delta/2$. Thus, any hyperplane inducing a different partition compared to $H(\v^m,b^m)$ has a margin at most $M-\delta/2$. Since $H(\v_m, b_m)$ is strictly feasible it therefore remains the unique maximum margin hyperplane in $F^\delta_{\mathrm{CL}}$. Observe now that for $0<h<\delta/2$ we have $H(\v^\star_{h,
\eta(h), L(h)}, b^\star_{h, \eta(h), L(h)}) \in F_{\mathrm{CL}}^\delta$, by Proposition \[prop:minim\]. In addition, by Theorem \[thm:convergence\], we know that the minimisers of $\hat I(\v, b)$ over $F^\delta_{\mathrm{CL}}$, say $H(\v_h^\delta, b_h^\delta)$, satisfy $$\lim_{h\to 0^+} \min \left\{\|(\v_h^\delta, b_h^\delta) - (\v^m, b^m)\|, \|(\v_h^\delta, b_h^\delta) + (\v^m, b^m)\|\right\} = 0.
$$ Now, since $H(\v^m, b^m)$ is strictly feasible $\exists \epsilon^\prime > 0$ s.t. $(\v, b) \in \mathbb{B}_{\epsilon^\prime}(\{(\v^m, b^m), -(\v^m, b^m)\})
\Rightarrow H(\v, b) \in F_{\mathrm{CL}}$. Then for any $0 < \epsilon <
\epsilon^\prime$ there exists $h^\prime > 0$ s.t. for $0<h<h^\prime$ both $(\v_h^\delta, b_h^\delta) \in \mathbb{B}_{\epsilon}(\{(\v^m, b^m), -(\v^m,
b^m)\}) \Rightarrow H(\v_h^\delta, b_h^\delta) \in F_{\mathrm{CL}}$ *and* $H(\v^\star_{h, \eta(h), L(h)}, b^\star_{h, \eta(h), L(h)}) \in
F_{\mathrm{CL}}^\delta$. Now for $H(\v, b) \in F^\delta_{\mathrm{CL}}\setminus
F_{\mathrm{CL}}$ we know that $\hat I(\v, b) < f_{\mathrm{CL}}(\v, b)$, whereas for $H(\v, b) \in F_{\mathrm{CL}}, \hat I(\v, b) = f_{\mathrm{CL}}(\v, b)$ and therefore the minimiser of $f_{\mathrm{CL}}(\v, b)$ must lie in the neighbourhood $\mathbb{B}_{\epsilon}(\{(\v^m, b^m), -(\v^m, b^m)\})$, and the result follows.

To illustrate the convergence of the MDH to the MMH we use the two-dimensional data set shown in Figure \[fig:MMHCon\]. The data is sampled from a mixture of two Gaussian distributions with equal covariance matrix. The MDH with respect to the true underlying density is $H\left( (1,-1),0
\right)$. A large margin separator is artificially introduced by removing a few observations in a narrow margin around a hyperplane different from $H\left(
(1,-1),0 \right)$. The margin is intentionally small to ensure that identifying the MMH is non-trivial. Figure \[fig:MMHCon\] illustrates the MDH solutions arising from the method for a decreasing sequence of bandwidths, $h$. Initially the MDH approximately coincides with the optimal MDH with respect to the true density of the Gaussian mixture. As $h$ decreases, the MDH approaches the MMH and for the smallest values of $h$ the two are indistinguishable.

![Convergence of the MDH to the maximum margin hyperplane for a decreasing sequence of bandwidth parameters, $h$.[]{data-label="fig:MMHCon"}](mmh_Gauss.pdf)

for Semi-Supervised Classification
----------------------------------

Denote the set of hyperplanes which correctly classify the labelled data by $F_{\mathrm{LB}}$. Under the assumption that $\exists H(\v, b) \in F_{\mathrm{LB}}\cap F_{\mathrm{CL}}$ with non-zero margin, we can show that, provided the parameter $\gamma$ does not shrink too quickly with $h$, the hyperplane that minimises $\fSL$ converges to the MMH contained in $F_{\mathrm{LB}} \cap
F_{\mathrm{CL}}$, where as before we assume that such an MMH is strictly feasible. To establish this result it is sufficient to show that there exists $h'>0$ such that for all $h \in (0,h^{\prime})$, the optimal hyperplane $H(\v^\star_{h,\eta,L,\gamma}, b^\star_{h, \eta,L,
\gamma})$ correctly classifies all the labelled examples. If this holds, then $\fSL(\v^\star_{h,\eta,L,\gamma}, b^\star_{h,\eta,L,\gamma}) = \fUL(\v^\star_{h,\eta,L,\gamma}, b^\star_{h,\eta,L,\gamma})$ for all sufficiently small $h$, and hence Lemma \[lem:ULConv\] can be applied to establish the result. The proof relies on the fact that the penalty terms associated with the known labels in Equation (\[eq:PenalisedSSL\]) are polynomials in $b$. Provided that $\gamma$ is bounded below by a polynomial in $h$, the value of the penalty terms for hyperplanes that do not correctly classify the labelled data dominate the value of the density integral as $h$ approaches zero. Therefore the optimal hyperplane must correctly classify the labelled data for small values of $h$.

Define $F_{\mathrm{LB}} = \{H(\v, b) \big \vert y_i(\v\cdot x_i - b) > 0, \forall
i=1, \dots, \ell\}$ and $F_{\mathrm{CL}} = \{H(\v, b) \big \vert \mu_{\v} -
\alpha\sigma_{\v} \leqslant b \leqslant \mu_{\v} + \alpha\sigma_{\v}\}$ and assume that $F_{\mathrm{SSC}} = F_{\mathrm{LB}} \cap F_{\mathrm{CL}} \not = \emptyset$ and that $\exists H(\v, b) \in 
F_{\mathrm{SSC}}$ with non-zero margin. For $h > 0$, let $L(h) = (e^{1/2}h^2\sqrt{2\pi})^{-1}$, $0 < \eta(h) \leqslant h$ and $\gamma(h) \geqslant h^r$ for some $r>0$. Then $\exists h^\prime > 0$ s.t. $h \in(0,h^\prime) \Rightarrow H(\v^\star_{h, \eta(h), L(h), \gamma(h)},
b^\star_{h, \eta(h), L(h), \gamma(h)} ) \in F_{\mathrm{LB}}$.

Consider $H(\v, b) \not \in F_{\mathrm{LB}}$. Then, $$\fSL(\v, b) \geqslant \frac{1}{n\sqrt{2\pi}h}\exp(-\nu_\star^2/2h^2) + \gamma(h)\nu_\star^{1+\epsilon} > \gamma(h)\nu_\star^{1+\epsilon},$$ where $\nu_\star > 0$ minimises $\frac{1}{n\sqrt{2\pi}h}\exp(-\nu^2/2h^2) + \gamma(h)\nu^{1+\epsilon}$. Therefore, $\nu_\star$ is the unique positive number satisfying, $$\begin{aligned}
\frac{1}{n\sqrt{2\pi}h}\exp\left(-\frac{\nu_\star^2}{2h^2}\right)\left(-\frac{\nu_\star}{h^2}\right) &+ (1+\epsilon)\gamma(h)\nu_\star^\epsilon = 0\\ &\Rightarrow \nu_\star^{1-\epsilon} = (1+\epsilon)\gamma(h)n\sqrt{2\pi}h^3\exp\left(\frac{\nu_\star^2}{2h^2}\right)\\
&\Rightarrow \nu_\star \geqslant \left((1+\epsilon)\gamma(h)n\sqrt{2\pi}h^3\right)^{1/1-\epsilon}.
$$ We therefore have, $$\begin{aligned}
\fSL(\v, b) &>& \gamma(h)\left((1+\epsilon)\gamma(h)n\sqrt{2\pi}h^3\right)^{\frac{1+\epsilon}{1-\epsilon}}\\
&=& K \gamma(h)^{\frac{2}{1-\epsilon}}h^{\frac{3(1+\epsilon)}{1-\epsilon}}\\
&\geqslant& K h^{\frac{2r+3(1+\epsilon)}{1-\epsilon}},
$$ where $K$ is a constant which can be chosen independent of $(\v, b)$. Finally, for any $H(\v^\prime, b^\prime) \in F_{\mathrm{SSC}}$ with non-zero margin, $\exists h^\prime > 0$ s.t. $$h \in(0,h^\prime) \Rightarrow \fSL(\v^\prime,
b^\prime) = \hat I(\v^\prime, b^\prime) < K
h^{\frac{2r+3(1+\epsilon)}{1-\epsilon}} < \fSL(\v, b).$$ Since $K$ is independent of $(\v, b)$, the result follows. The final set of inequalities holds since the hyperplane $H(\v^\prime, b^\prime)$ is assumed to have non-zero margin, say $M_{\v^\prime, b^\prime} > 0$, and hence $\hat I(\v^\prime, b^\prime) \leqslant \frac{1}{h\sqrt{2\pi}}\exp\{-M_{\v^\prime, b^\prime}/2h^2\}$, which tends to zero faster than any polynomial in $h$.

Estimation of Minimum Density Hyperplanes {#sec:Methodology}
=========================================

In this section we discuss the computation of MDHs. We first investigate the continuity and differentiability properties required to optimise the projection indices $\phiUL(\v)$ and $\phiSSL(\v)$.

Since the domain of both projection indices, $\phiUL(\v)$ and $\phiSSL(\v)$, is the boundary of the unit-sphere in $\R^d$ it is more convenient to express $\v$ in terms of spherical coordinates, $$\label{eq:Spherical}
v_i(\theta)  = \left\{\begin{array}{ll}\cos(\theta_i) \prod_{j=1}^{i-1}
\sin(\theta_j),& i=1,\ldots,d-1  \\
\prod_{j=1}^{d-1} \sin(\theta_j), & i = d, \end{array}\right.
$$ where $\theta \in \Th = [0, \pi]^{d-2} \times [0,2\pi]$ is called the [*projection angle*]{}. Using spherical coordinates renders the domain, $\Th$, convex and compact, and reduces dimensionality by one.

As the following discussion applies to both $\phiUL(\v)$ and $\phiSSL(\v)$ we denote a generic projection index $\phi:\Theta \to \R$, and the associated set of minimisers, as, $$\begin{aligned}
\label{OpVal}
\phi(\theta) &= \min_{b \in A} f(\v(\theta),b),\\
\hY(\theta) &= \left\{ b \in A\, \big| \, f(\v(\theta), b) = \phi(\theta) \right\},
$$ where $f(\v(\theta), b)$ is continuously differentiable, $A \subset \R$ is compact and convex, and the correspondence $\hY(\theta)$ gives the set of global minimisers of $f(\v(\theta),b)$ for each $\theta$. The definition of $A$ is not critical in our formulation. Setting, $$A \supset \left[\min_{\v \in \B}\{\mu_\v\}-\alpha \sigma_{\textrm{pc}_1} -\eta, \max_{\v \in \B}\{\mu_\v\}+\alpha \sigma_{\textrm{pc}_1} + \eta \right],
$$ where $\sigma^2_{\textrm{pc}_1}$ is the variance of the projections along the first principal component, ensures that the set of hyperplanes that satisfy the constraint of Equation (\[constr1\]) will be a subset of $A$ for all $\v$.

Berge’s maximum theorem [@Berge1963; @Polak1987], establishes the continuity of $\phi(\theta)$ and the upper-semicontinuity (u.s.c.) of the correspondence $\hY(\theta)$. Theorem 3.1 in [@Polak1987] enables us to establish that $\phi(\theta)$ is locally Lipschitz continuous. Using Theorem 4.13 of @BonnansS2000 we can further show that $\phi(\theta)$ is directionally differentiable everywhere. The directional derivative at $\theta$ in the direction $\nu$ is given by, $$\label{OptDif}
d \phi(\theta; \nu) = \min_{b \in \hY(\theta)} D_{\theta} f(\v(\theta), b)
\cdot \nu,
$$ where $D_{\theta}$ denotes the derivative with respect to $\theta$. It is clear from Equation (\[OptDif\]) that $\phi(\theta)$ is differentiable if $D_{\theta} f(\v(\theta), b)$ is the same for all $b \in \hY(\theta)$. If $B(\theta)$ is a singleton then this condition is trivially satisfied and $\phi(\theta)$ is continuously differentiable at $\theta$.

It is possible to construct examples in which $B(\theta)$ is not a singleton. However, with the exception of contrived examples, our experience with real and simulated data sets indicates that when $h$ is set through standard bandwidth selection rules $B(\theta)$ is almost always a singleton over the optimisation path.

Suppose $B(\theta)$ is a singleton for almost all $\theta \in \Th$. Then $\phi(\theta)$ is continuously differentiable almost everywhere.

The result follows immediately from the fact that if $B(\theta) = \{b\}$ is a singleton, then the derivative $D\phi(\theta) = D_\theta f(\v(\theta), b)$, which is continuous.

@Wolfe1972 has provided early examples of how standard gradient-based methods can fail to converge to a local optimum when used to minimise nonsmooth functions. In the last decade a new class of nonsmooth optimisation algorithms has been developed based on gradient sampling [@BurkeLO2006]. Gradient sampling methods use generalised gradient descent to find local minima. At each iteration points are randomly sampled in a radius $\varepsilon$ of the current candidate solution, and the gradient at each point is computed. The convex hull of these gradients serves as an approximation of the $\varepsilon$-Clarke generalised gradient [@BurkeLO2002a]. The minimum element in the convex hull of these gradients is a descent direction. The gradient sampling algorithm progressively reduces the sampling radius so that the convex hull approximates the Clarke generalised gradient. When the origin is contained in the Clarke generalised gradient there is no direction of descent, and hence the current candidate solution is a local minimum. Gradient sampling achieves almost sure global convergence for functions that are locally Lipschitz continuous and almost everywhere continuously differentiable. It is also well documented that it is an effective optimisation method for functions that are only locally Lipschitz continuous.

Computational Complexity
------------------------

In this subsection we analyse the computational complexity of MDP$^2$. At each iteration the algorithm projects the data sample onto $\v(\theta)$ which involves $\O(n d)$ operations. To compute the projection index, $\phi(\theta)$, we need to minimise the penalised density integral, $f(\v(\theta),b)$. This can be achieved by first evaluating $f(\v(\theta),b)$ on a grid of $m$ points, to bracket the location of the minimiser, and then applying bisection to compute the minimiser(s) within the desired accuracy. The main computational cost of this procedure is due to the first step which involves $m$ evaluations of a kernel density estimator with $n$ kernels. Using the improved fast Gauss transform [@Morariu08] this can be performed in $\O(m+n)$ operations, instead of $\O(mn)$. Bisection requires $\O(-\log_2
\varepsilon)$ iterations to locate the minimiser with accuracy $\varepsilon$.

If the minimiser of the penalised density integral $b^\star = \operatorname*{arg\,min}_{b \in A}
f(\v(\theta),b)$, is unique the projection index is continuously differentiable at $\theta$. To obtain the derivative of the projection index it is convenient to define the projection function, $P(\v) = \left( \x_1 \cdot \v, \ldots, \x_n \cdot \v \right)^\top.$ An application of the chain rule yields, $$\label{Der}
d_\theta \phi = D_\theta f(\v(\theta), b^\star) = D_P f(\v(\theta), b^\star)
D_\v P D_\theta \v
$$ where the derivative of the projections of the data sample with respect to $\v$ is equal to the data matrix, $D_\v P = (\x_1, \ldots, \x_n)^{\top}$; and $D_\theta \v$ is the derivative of $\v$ with respect to the projection angle, which yields a $d\times (d-1)$ matrix. The computation of the derivative therefore requires $\O(d(n+d))$ operations.

The original GS algorithm requires $\O(d)$ gradient evaluations at each iteration which is costly. @CurtisQ2013 have developed an adaptive gradient sampling algorithm that requires $\O(1)$ gradient evaluations in each iteration. More recently, @LewisO2013 have strongly advocated that for the minimisation of nonsmooth, nonconvex, locally Lipschitz functions, a simple BFGS method using inexact line searches is much more efficient in practice than gradient sampling, although no convergence guarantees have been established for this method. BFGS requires a single gradient evaluation at each iteration and a matrix vector operation to update the Hessian matrix approximation. In our experiments we use the BFGS algorithm.

Experimental Results {#sec:Exper}
====================

In this section we assess the empirical performance of MDHs for clustering and semi-supervised classification. We compare performance with existing state-of-the-art methods for both problems on the following 14 benchmark data sets: Banknote authentication (banknote), Breast Cancer Wisconsin original (br. cancer), Forest type mapping (forest), Ionosphere, Optical recognition of handwritten digits (optidigits), Pen-based recognition of hand-written digits (pendigits), Seeds, Smartphone-Based Recognition of Human Activities and Postural Transitions (smartphone), Statlog Image Segmentation (image seg.), Statlog Landsat Satellite (satellite), Synthetic control chart time series (synth control), Congressional voting records (voting), Wine, and Yeast cell cycle analysis (yeast). Details of these data sets, in terms of their size, $n$, dimensionality, $d$ and number of clusters, $c$, can be seen in Table \[tb:datasets\].

                        $n$   $d$   $c$
  ----------------- ------- ----- -----
       banknote$^a$    1372     4     2
     br. cancer$^a$     699     9     2
         forest$^a$     523    27     4
     ionosphere$^a$     351    33     2
     optidigits$^a$    5618    64    10
      pendigits$^a$   10992    16    10
          seeds$^a$     210     7     3
    smartphone $^a$   10929   561    12
     image seg.$^a$    2309    18     7
      satellite$^a$    6435    36     6
          synth$^a$     600    60     6
         voting$^a$     435    16     2
           wine$^a$     178    13     3
          yeast$^b$     698    72     5

  : Details of benchmark data sets: size ($n$), dimensionality ($d$), number of clusters ($c$).[]{data-label="tb:datasets"}

[a. UCI machine learning repository <https://archive.ics.uci.edu/ml/datasets.html>\
b. Stanford Yeast Cell Cycle Analysis Project <http://genome-www.stanford.edu/cellcycle/> ]{}

Clustering {#ssec:expUL}
----------

Since an MDH yields a bi-partition of a data set rather than a complete clustering, we propose two measures to assess the quality of a binary partition of a data set containing an arbitrary number of clusters. Both take values in $[0, 1]$ with larger values indicating a better partition. These measures are motivated by the fact that a good binary partition should (a) avoid dividing clusters between elements of the partition, and (b) be able to discriminate at least one cluster from the rest of the data. To capture this we modify the cluster labels of the data by assigning each cluster to the element of the binary partition which contains the majority of its members. In the case of a tie the cluster is assigned to the smaller of the two partitions. We thus merge the true clusters into two aggregate clusters, $C_{1}$ and $C_{2}$.

The first measure we use is the binary V-measure which is simply the V-measure [@Rosenberg2007] computed on $C_{1}, C_{2}$ with respect to the binary partition, which we denote $\Pi_1, \Pi_2$. The V-measure is the harmonic mean of homogeneity and completeness. For a data set containing clusters $C_1,\ldots,C_c$, partitioned as $\Pi_1,\ldots,\Pi_k$, homogeneity is defined as the conditional entropy of the cluster distribution within each partition, $\Pi_i$. Completeness is symmetric to homogeneity and measures the conditional entropy of each partition within each cluster, $C_j$. An important characteristic of the V-measure for evaluating binary partitions is that if the distribution of clusters within each partition is equal to the overall cluster distribution in the data set then the V-measure is equal to zero [@Rosenberg2007]. This means that if an algorithm fails to distinguish the majority of any of the clusters from the remainder of the data, the binary V-measure returns zero performance. Other evaluation metrics for clustering, such as purity and the Rand index, can assign a high value to such partitions.

To define the second performance measure we first determine the number of correctly and incorrectly classified samples. The error of a binary partition, $\mbox{E}(\Pi_1, \Pi_2)$, given in Equation (\[eq:error\]), is defined as the number of elements of each aggregate cluster which are not in the same partition as the majority of their original clusters. In contrast, the success of a partition, $\mbox{S}(\Pi_1, \Pi_2)$, Equation (\[eq:success\]), measures the number of samples which are in the same partition as the majority of their original clusters. The Success Ratio, $\mbox{SR}(\Pi_1, \Pi_2)$, Equation (\[eq:SR\]), captures the extent to which the majority of at least one cluster is well-distinguished from the rest of the data. $$\begin{aligned}
\mbox{E}(\Pi_1, \Pi_2) & = \min\left\{\vert \Pi_1 \cap C_1 \vert + \vert \Pi_2
    \cap C_2 \vert, \vert \Pi_1 \cap C_2 \vert + \vert \Pi_2 \cap C_1 \vert\right\}, \label{eq:error}\\
\mbox{S}(\Pi_1, \Pi_2) & = \min\left\{\max\left\{\vert\Pi_1 \cap C_1 \vert, \vert
\Pi_1 \cap C_2 \vert \right\}, \max\left\{\vert \Pi_2 \cap C_1 \vert, \vert \Pi_2 \cap C_2 \vert \right\}\right\}, \label{eq:success}\\
\mbox{SR}(\Pi_1, \Pi_2) & = \frac{ \mbox{S}(\Pi_1, \Pi_2) }{ \mbox{S}(\Pi_1, \Pi_2) + \mbox{E}(\Pi_1, \Pi_2)}. \label{eq:SR}
$$ The Success Ratio takes the value zero if an algorithm fails to distinguish the majority of any cluster from the remainder of the data.

### Parameter Settings for 

The two most important settings for the performance of the proposed approach are the initial projection direction, and the choice of $\alpha$, which controls the width of the interval $F(\v)$ within which the optimal hyperplane falls. Despite the ability of the formulation to mitigate the effect of local minima of the projected density, the problem remains non-convex and local minima in the projection index can still lead to suboptimal performance. We have found that this effect is amplified in general when either or both the number of dimensions, and the number of high density clusters in the data set is large. To better handle the effect of local optima, we use multiple initialisations and select the MDH that maximises the *relative depth* criterion, defined in Equation (\[eq:relDep\]). The relative depth of an MDH, $H(\v, b)$, is defined as the smaller of the relative differences in the density on the MDH and its two adjacent modes in the projected density, $$\label{eq:relDep}
\mbox{RelativeDepth}(\v, b) = \frac{ \min\left\{ \hat{I}(\v, m_l), \hat{I}(\v, m_r) \right\} - \hat{I}(\v, b)}{\hat{I}(\v, b)}
$$ where $m_l$ and $m_r$ are the two adjacent modes in the projected density on $\v$. If an MDH does not separate the modes of the projected density, then its relative depth is set to zero, signalling a failure of to identify a meaningful bi-partition. The relative depth is appealing because it captures the fact that a high quality separating hyperplane should have a low density integral, and separate well the modes of the projected density. Note also that the relative depth is equivalent to the inverse of a measure used to define cluster overlap in the context of Gaussian mixtures [@Aitnouri2000]. In all the reported experiments we initialise to the first and second principal component and select the MDH with the largest relative depth. For the data sets listed above it was never the case that both initialisations led to MDHs with zero relative depth.

The choice of $\alpha$ determines the trade-off between a balanced bi-partition and the ability to discover lower density hyperplanes. The difficulties associated with choosing this parameter are illustrated in Figure \[fig:alpha\]. In each sub-figure the horizontal axis is the candidate projection vector, $\v$, while the right vertical axis is the direction of maximum variability orthogonal to $\v$. Points correspond to projections of the data sample onto this two-dimensional space, while colour indicates cluster membership. The solid line depicts the projected density on $\v$, while the dotted line depicts the penalised function, $\fUL(\v, \cdot)$. The scale of both functions is depicted on the left vertical axis. The solid vertical line indicates the MDH along $\v$. Setting $\alpha$ to a large value can cause to focus on hyperplanes that have low density because they partition only a small subset of the data set as shown in Figure \[fig:alpha1\]. In contrast smaller values of $\alpha$ may cause the algorithm to disregard valid lower density hyperplane separators (see Figure \[fig:alpha2\]), or for the separating hyperplane to not be a local minimiser of the projected density (see Figure \[fig:alpha3\]).

Rather than selecting a single value for $\alpha$ we recommend solving repeatedly for an increasing sequence of values in the range $\{\alpha_{\min}, \alpha_{\max}\}$, where each implementation beyond the first is initialised using the solution to the previous. Setting $\alpha_{\min}$ close to zero forces to seek low density hyperplanes that induce a balanced data partition. This tends to find projections which display strong multimodal structure, yet prevents convergence to hyperplanes that have low density because they partition a few observations, as in the case shown in Figure \[fig:alpha1\]. Increasing $\alpha$ progressively fine-tunes the location of the MDH. To avoid sensitivity to the value of $\alpha_{\max}$ (set to 0.9) the output of the algorithm is the last hyperplane that corresponds to a minimiser of the projected density. Figure \[fig:opti\] illustrates this approach using the optical recognition of handwritten digits data set from the UCI machine learning repository [@uci]. Figure \[fig:opt1\] depicts the projected density on the initial projection direction, which in this case is the second principal component. As shown, the density is unimodal and the clusters are not well separated along this vector. Although not shown, if a large value of $\alpha$ is used from the outset, will identify a vector along which the projected density is unimodal and skewed. Figure \[fig:opt2\] shows that after five iterations with $\alpha=10^{-2}$ has identified a projection vector with bimodal density. In subsequent iterations the two modes become more clearly separated, Figure \[fig:opt3\], while increasing $\alpha$ enables to locate an MDH that corresponds to a minimiser of $\hat{I}(\v, b)$, as illustrated in Figure \[fig:opt4\].

\

In all experiments we set the bandwidth parameter to $h = 0.9
\hat{\sigma}_{\textrm{pc}_1}n^{-1/5}$, where $\hat{\sigma}_{\textrm{pc}_1}$ is the estimated standard deviation of the data projected onto the first principal component. This bandwidth selection rule is recommended when the density being approximated is assumed to be multimodal [@Silverman1986]. The parameter $\eta$ controls the distance between the minimisers of $\operatorname*{arg\,min}_{b \in \R} \fUL(\v,b)$ and $\operatorname*{arg\,min}_{b \in
F(\v)} \hat{I}(\v,b)$, while larger values of $\epsilon$ increase the smoothness of the penalised function $\fUL$. Values of $\eta$ close to zero affect the numerical stability of the one-dimensional optimisation problem, due to the term $\frac{L}{\eta^\epsilon}$ in $\fUL$ becoming very large. We used $\eta=10^{-2}$ and $\epsilon=1-10^{-6}$ to avoid numerical instability. Beyond these numerical problems the values of $\eta$ and $\epsilon$ do not affect the solutions obtained through MDP$^2$.

### Performance Evaluation

We compare the performance of for clustering with the following methods:

1.  $k$-means++ [@ArthurV2007], a version of $k$-means that is guaranteed to be $\mathcal{O}(\log k)$-competitive to the optimal $k$-means clustering.

2.  The adaptive linear discriminant analysis guided $k$-means (LDA-$k$m) [@DingL2007]. LDA-$k$m attempts to discover the most discriminative linear subspace for clustering by iteratively using $k$-means, to assign labels to observations, and LDA to identify the most discriminative subspace.

3.  The principal direction divisive partitioning (PDDP) [@Boley1998], and the density-enhanced PDDP (dePDDP) [@TasoulisTP2010]. Both methods project the data onto the first principal component. PDDP splits at the mean of the projections, while dePDDP splits at the lowest local minimum of the one-dimensional density estimator.

4.  The iterative support vector regression algorithm for MMC [@ZhangTK2009] using the inner product and Gaussian kernel, iSVR-L and iSVR-G respectively. Both are initialised with the output of 2-means++.

5.  Normalised cut spectral clustering (SCn) [@NgJW2002] using the Gaussian affinity function, and the automatic bandwidth selection method of [@Zelnik2004]. This choice of kernel and bandwidth produced substantially better performance than alternative choices considered. For data sets that are too large for the eigen decomposition of the Gram matrix to be feasible we employed the Nystr[ö]{}m method [@FowlkesBSM2004].

We also considered the density-based clustering algorithm PdfCluster [@MenardiA2014], but this algorithm could not be executed on the larger data sets and so its performance is not reported in this paper. With the exception of SCn and iSVR-G, the methods considered bi-partition the data through a hyperplane in the original feature space. For the 2-means and LDA-2m algorithm the hyperplane separator bisects the line segment joining the two centroids. iSVR-L directly seeks the maximum margin hyperplane in the original space, while iSVR-G seeks the maximum margin hyperplane in the feature space defined by the Gaussian kernel. PDDP and dePDDP use a hyperplane whose normal vector is the first principal component. PDDP uses a fixed split point while dePDDP uses the hyperplane with minimum density along the fixed projection direction.

Table \[tbl:cluster\] reports the performance of the considered methods with respect to the success ratio (SR) and the binary V-measure (V-m) on the fourteen data sets. In addition Figures \[fig:clMetr\] and \[fig:clReg\] provide summaries of the overall performance on all data sets using boxplots of the raw performance measures as well as the associated *regret*. The regret of an algorithm on a given data set is defined as the difference between the best performance attained on this data set and the performance of this algorithm. By comparing against the best performing clustering algorithm regret accommodates for differences in difficulty between clustering problems, while also making use of the magnitude of performance differences between algorithms. The distribution of performance with respect to both SR and V-m is negatively skewed for most methods, and as a result the median is higher than the mean (indicated with a red dot).

It is clear from Table \[tbl:cluster\] that no single method is consistently superior to all others, although achieves the highest or tied highest performance on seven data sets (more than any other method). More importantly is among the best performing methods in almost all cases. This fact is better captured by the regret distributions in Figure \[fig:clReg\]. Here we see that the average, median, and maximum regret of is substantially lower than any of the competing methods. In addition achieves the highest mean and median performance with respect to both SR and V-m, while also having much lower variability in performance when compared with most other methods.

Pairwise comparisons between and other methods reveal some less obvious facts. SCn achieves higher performance than in more examples (six) than any other competing method, however it is much less consistent in its performance, obtaining very poor performance on five of the data sets. The iSVR maximum margin clustering approach is arguably the closest competitor to MDP$^2$. iSVR-L and iSVR-G achieve the second and third highest average performance with respect to V-m and SR respectively. The PDDP algorithm is the second best performing method on average with respect to SR, but performs poorly with respect to V-m. The density enhanced variant, dePDDP, performs on average much worse than MDP$^2$. This approach is similarly motivated by obtaining hyperplanes with low density integral, and its low average performance indicates the usefulness of searching for high quality projections as opposed to always using the first principal component. Finally, neither of the $k$-means variants appears to be competitive with in general.

Semi-Supervised Classification {#ssec:expSSL}
------------------------------

In this section we evaluate MDHs for semi-supervised classification. We compare MDHs against three state-of-the-art semi-supervised classification methods: Laplacian Regularised Support Vector Machines (LapSVM) [@BelkinNS2006], Simple Semi-Supervised Learning (SSSL) [@JiYLJH2012], and Correlated Nystr[" o]{}m Views (XNV) [@mcwilliams2013]. For all methods the inner product kernel was used to render the resulting classifiers linear, and thereby comparable to our method. As the MDH is asymptotically equivalent to a linear §3VM we also considered the continuous formulation for the estimation of a §3VM proposed by [@ChapelleZ2005]. These results are omitted as this method was not competitive on any of the considered data sets.

### Parameter Settings for 

The existence of a few labelled examples enables an informed initialisation of MDP$^2$. We consider the first and second principal components as well as the weight vector of a linear SVM trained on the labelled examples only, and initialise with the vector that minimises the value of the projection index,$~\phiSSL$. The penalty parameter $\gamma$ is first set to $0.1$ and with this setting $\alpha$ is progressively increased in the same way as for clustering. After this, $\alpha$ is kept at $\alpha_{\max}$ and $\gamma$ is increased to 1 and then 10. Thus the emphasis is initially on finding a low density hyperplane with respect to the marginal density $\hat{p}(\x)$. As the algorithm progresses the emphasis on correctly classifying the labelled examples increases, so as to obtain a hyperplane with low training error within the region of low density already determined.

### Performance Evaluation

To assess the effect on performance of the number of labelled examples, $\ell$, we consider a range of values. We compare the methods using the subset of data sets used in the previous section in which the size of the smallest class exceeds 100. In total eight data sets are used. For each value of $\ell$, 30 random partitions into labelled and unlabelled data are considered. As classes are balanced in the data sets considered, performance is measured only in terms of classification error on the unlabelled data. For data sets with more than two classes all pairwise combinations of classes are considered and aggregate performance is reported.

median (—[$\triangle$]{}—), LapSVM median (), SSSL median (), XNV median ([—[$+$]{}—]{}), with corresponding interquartile ranges given by shaded regions.

Figure \[ssl:regret1\] provides plots of the median and interquartile range of the classification error for values of $\ell$ between 5 and 100 for the four data sets with two classes. Overall appears to be most competitive when the number of labelled examples is small. In addition, is comparable with the best performing method in almost every case. The only exception is the ionosphere data set where LapSVM outperforms for all values of $\ell$. Figure \[ssl:regret2\] provides plots of the median and interquartile range of the aggregate classification error on data sets containing more than two classes. As these data sets are larger we consider up to 300 labelled examples. Note that the interquartile range for XNV is not depicted for the satellite data set. The variability of performance of XNV on this data set was so high that including the interquartile range would obscure all other information in the figure. exhibits the best performance overall, and obtains the lowest median classification error, or tied lowest, for all data sets and values of $\ell$.

median (—[$\triangle$]{}—), LapSVM median (), SSSL median (), XNV median ([—[$+$]{}—]{}), with corresponding interquartile ranges given by shaded regions.

Summary of Experimental Results
-------------------------------

We evaluated the performance of the formulation for finding MDHs for both clustering and semi-supervised classification, on a large collection of benchmark data sets, and in comparison with state-of-the-art methods for both problems.

For clustering, we found that no single method was consistently superior to all others. This is a result of the vastly differing nature of the data sets in terms of size, dimensionality, number and shape of clusters, etc. achieved the best performance on more data sets than any of the competing methods, and importantly was competitive with the best performing method in almost every data set considered. All other methods performed poorly in at least as many examples. Boxplots of both the raw performance and performance regret, which measures the difference between each method and the best performing method on each data set, allowed us to summarise the comparative performance of the different methods across data sets. The mean and median raw performance of is substantially higher than the next best performing method, and the regret is also substantially lower.

In the case of semi-supervised classification it was apparent that is extremely competitive when the number of labelled examples is (very) small, but that in some cases its performance does not improve as much as that of the other methods considered, when the labelled examples become more abundant. Our experiments suggest that overall is very competitive with the state-of-the-art for semi-supervised classification problems.

Conclusions {#sec:concl}
===========

We proposed a new hyperplane classifier for clustering and semi-supervised classification. The proposed approach is motivated by determining low density linear separators of the high-density clusters within a data set. This is achieved by minimising the integral of the empirical density along the hyperplane, which is computed through kernel density estimation. To the best of our knowledge this is the first direct implementation of the low density separation assumption that underlies high-density clustering and numerous influential semi-supervised classification methods. We show that the minimum density hyperplane is asymptotically connected with maximum margin hyperplane, thereby establishing an important link between the proposed approach, maximum margin clustering, and semi-supervised support vector machines.

The proposed formulation allows us to evaluate the integral of the density on a hyperplane by projecting the data onto the vector normal to the hyperplane, and estimating a univariate kernel density estimator. This enables us to apply our method effectively and efficiently on data sets of much higher dimensionality than is generally possible for density based clustering methods. To mitigate the problem of convergence to locally optimal solutions we proposed a projection pursuit formulation.

We evaluated the minimum density hyperplane approach on a large collection of benchmark data sets. The experimental results obtained indicate that the method is competitive with state-of-the-art methods for clustering and semi-supervised classification. Importantly the performance of the proposed approach displays low variability across a variety of data sets, and is robust to differences in data size, dimensionality, and number of clusters. In the context of semi-supervised classification, the proposed approach shows especially good performance when the number of labelled data is small.

Proof of Theorem \[thm:convergence\] {#app:thm:convergence}
====================================

Before proving Theorem \[thm:convergence\] we require the following two technical lemmata which establish some algebraic properties of the maximum margin hyperplane. The following lemma shows that any hyperplane orthogonal to the maximum margin hyperplane results in a different partition of the support points of the maximum margin hyperplane. The proof relies on the fact that if this statement does not hold then a hyperplane with larger margin exists which is a contradiction. Figure \[fig:MMHorthog\] provides an illustration of why this result holds. (a) Any hyperplane orthogonal to MMH generates a different partition of the support points of MMH, e.g., the point highlighted in red in (b) is grouped with the lower three by the dotted line but with the upper two by the solid line, the MMH. If an orthogonal hyperplane *can* generate the same partition (c), then a larger margin hyperplane than the proposed MMH exists (d).

\
\
Proposed MMH [**—–**]{}, Orthogonal hyperplane [**- - -**]{}, Hyperplane with larger\
margin , Regular points $\bigcirc$, Support points $\pmb{\bigcirc}$,\
Differently assigned support point

\[lem:MMHorthog\]

Suppose there is a unique hyperplane in $F$ with maximum margin, which can be parameterised by $(\v^m, b^m) \in \B \times \R$. Let $M = \margin\, H(\v^m, b^m)$, $C^+ = \{\x \in \X \,\vert \, \v^m \cdot \x - b^m = M\}$ and $C^- = \{\x \in \X \,
\vert \, b^m - \v^m \cdot \x = M\}$. Then, $\forall \w \in \mathrm{Null}(\v^m)$, $c \in \R$ either $\min\{\w\cdot \x - c
\, \vert \, \x \in C^+\}\leqslant 0$, or $\max\{\w\cdot \x - c \, \vert \, \x \in
C^-\}\geqslant 0$.

Suppose the result does not hold, then $\exists (\w, c)$ with $\|\w\| = 1, \w
\cdot \v^m = 0$ and $\min\{\w \cdot \x - c \vert \x \in C^+\} > 0$ and $\max\{\w \cdot
\x - c \vert \x \in C^-\} < 0$. Let $m = \min\{ |\w\cdot \x - c | \,\big \vert\, \x \in C^+ \cup C^-\}$. Define $\lambda = \frac{m}{2M} < 1$. Define $\u = \frac{1}{\sqrt{\lambda^2 + (1-\lambda)^2}}(\lambda \w + (1-\lambda)
\v^m)$ and $d = \frac{\lambda c + (1-\lambda) b^m}{\sqrt{\lambda^2 +
(1-\lambda)^2}}$. By construction $\|\u\| = 1$. For any $\x_+ \in C^+$ we have, $$\begin{aligned}
\u \cdot \x_+ - d &= \frac{\lambda(\w \cdot \x_+ - c) + (1-\lambda)(\v^m \cdot \x_+ -
b^m)}{\sqrt{\lambda^2 + (1-\lambda)^2}}\\
&\geqslant \frac{\lambda m + (1-\lambda)M}{\sqrt{\lambda^2 + (1-\lambda)^2}}\\
&= \frac{m^2 + 2M^2 - Mm}{\sqrt{m^2 + (2M-m)^2}}\\
&>M.
$$ Similarly one can show that $d - \u \cdot \x_- > M$ for any $\x_- \in C^-$, meaning that $(\u,d)$ achieves a larger margin on $C^+$ and $C^-$ than $(\v^m, b^m)$, a contradiction.

The next lemma uses the above result to provide an upper bound on the distance between pairs of support points projected onto any vector, in terms of the angle between that vector and $\v^m$.

\[lm:2Mv\]

Suppose there is a unique hyperplane in $F$ with maximum margin, which can be parameterised by $(\v^m, b^m) \in \B \times \R$. Define $M =
\margin\, H(\v^m,b^m)$, $C^+ = \{\x \in \X \vert \v^m \cdot \x - b^m = M\}$, and $C^- = \{\x \in
\X \vert  b^m-\v^m\cdot \x = M\}$. There is no vector $\w \in \R^d$ for which $\w\cdot \x_+ - \w\cdot \x_- > 2M \v^m \cdot \w$ for all pairs $\x_+ \in C^+, \x_- \in C^-$.

Suppose such a vector exists. Define $\w' = \w - (\v^m\cdot \w)\v^m$. By construction $\w' \in
\mathrm{Null}(\v^m)$. For any pair $\x_+ \in C^+, \x_- \in C^-$ we have $$\begin{aligned}
\w' \cdot x_+ - \w' \cdot x_- &= \w \cdot x_+ - (\v^m \cdot \w) \v^m \cdot x_+ - \w
    \cdot x_- +  (\v^m \cdot \w) \v^m \cdot x_-\\
&> \v^m\cdot \w(2M - \v^m\cdot x_+ + b^m - b^m + \v^m \cdot x_-)\\
&= 0.  \end{aligned}$$ Define $c: = \frac{1}{2}(\min \{ \w' \cdot \x_ + \big \vert \x_+ \in C^+\} + \max
\{ \w' \cdot \x_- \big \vert \x_- \in C^-\})$. Then $\min\{\w'\cdot \x_+ - c \vert
\x_+ \in C^+\} > 0$ and $\max\{\w' \cdot \x_- - c \vert \x_- \in C^-\} < 0$, a contradiction.

We are now in a position to provide the main proof of this appendix. The theorem states that if the maximum margin hyperplane is unique, and can be parameterised by $(\v^m, b^m) \in \B\times \R$, then $$\lim_{h \to 0^+} \min \left\{\|(\v^\star_h, b^\star_h) - (\v^m, b^m)\|, \|(\v^\star_h, b^\star_h) + (\v^m, b^m)\|\right\} = 0,$$ where $\left\{H(v^\star_h, b^\star_h)\right\}_h$ is any collection of minimum density hyperplanes indexed by their bandwidth $h > 0$.

[[**of Theorem \[thm:convergence\]**]{}]{}

Define $M = \margin\, H(\v^m, b^m)$, $C^+ = \{\x \in \X \,\vert \, \v^m \cdot \x - b^m = M\}$ and $C^- = \{\x \in \X \, \vert \, b^m - \v^m \cdot \x = M\}$. Let $B = \max\{\| \x\| \, \big \vert \, \x \in \X\}$. Take any $\epsilon > 0$ and set $0<\delta$ to satisfy $\frac{2 \delta}{M}(1+B^2) + 2B\delta^{3/2}\sqrt{\frac{2}{M}} + \delta^2 = \epsilon^2$. Now, suppose $(\w, c) \in \B \times \R$ satisfies, $$\w \cdot \x_+ - c > M- \delta, \ \forall \x_+ \in C^+ \mbox{ and } c - \w \cdot \x_- > M - \delta, \ \forall \x_- \in C^-. 
$$ By Lemma \[lm:2Mv\] we know that $\exists \x_+ \in C^+, \x_- \in C^-$ s.t. $\w \cdot \x_+ -  \w\cdot \x_- \leqslant 2M \v^m \cdot \w$. Thus $$\begin{aligned}
\v^m \cdot \w &\geqslant \frac{\w\cdot \x_+ - \w \cdot \x_-}{2M}\\
&= \frac{\w\cdot \x_+ -c + c- \w \cdot \x_-}{2M}\\
&> \frac{2M - 2 \delta}{2M} = 1 - \frac{\delta}{M}.
$$ Thus $\|\v^m - \w \|^2 < \frac{2\delta}{M}$. Now, for each $\x_+ \in C^+, \v^m \cdot \x_+ - b = M$ and for each $\x_- \in C^-, b - \v^m \cdot \x_- = M$. Thus for any such $\x_+, \x_-$ we have, $$\begin{aligned}
&M - \delta + \w \cdot \x_- < c < \w \cdot \x_+ - M + \delta,\\
& b^m - \v^m \cdot \x_- - \delta + \w \cdot \x_- < c < \w \cdot \x_+ - \v^m \cdot \x_+ + b^m + \delta,\\
& b^m-\delta - (\v^m - \w)\cdot \x_- < c < b^m + \delta + (\w - \v^m)\cdot \x_+,\\
& b^m-\delta - B\|\v^m - \w\| < c < b^m+\delta + B\|\w-\v^m\|, \\
& | c - b^m | < \left| \delta + B\|\w-\v^m\| \right|.
$$ We can now bound the distance between $(\w,c)$ and $(\v^m,b^m)$, $$\begin{aligned}
\|(\v^m, b^m) - (\w, c)\|^2 &=& \|\v^m-\w\|^2 + \vert b^m -c\vert ^2\\
&<& \|\v^m-\w\|^2(1+B^2) + 2B\delta\|\v^m-\w\| + \delta^2\\
&<& \frac{2 \delta}{M}(1+B^2) + 2B\delta\sqrt{\frac{2 \delta}{M}} + \delta^2\\
&=& \epsilon^2.
$$

We have shown that for any hyperplane $H(\w, c)$ that achieves a margin larger than $M-\delta$ on the support points of the maximum margin hyperplane, $\x \in C^+ \cup C^-$, the distance between $(\w,c)$ and $(\v^m,b^m)$ is less than $\epsilon$. Equivalently, any hyperplane $H(\w,
c)$ such that $\| (\w,c) - (\v^m, b^m)\| > \epsilon$ has a margin less than $M-\delta$, as $\min\left\{ | \w \cdot \x - c| \, \big\vert \, \x \in C^+ \cup C_- \right\} < M-\delta$. By symmetry, the same holds for any $(\w,c)$ within distance $\epsilon$ of $(-\v^m, -b^m)$.

By Lemma \[lm:samedivision\] $\exists h_1>0$ such that for all $h \in (0,h_1)$, the minimum density hyperplane for $h$, $H(\v^\star_h, b^\star_h)$, induces the same partition of $\X$ as the maximum margin hyperplane, $H(\v^m, b^m)$. By Lemma \[lm:bddmargin\] $\exists h_2 > 0$ such that for all $h \in(0,h_2)$, $\margin\, H(\v^{\star}_h, b^{\star}_h) > M-\delta$. Therefore for $h \in (0,\min\{h_1, h_2\})$, $\min\left\{ \|(\v^{\star}_h, b^{\star}_h) - (\v^m, b^m)\|,
\| (\v^{\star}_h, b^{\star}_h) +(\v^m, b^m)\| \right\} < \epsilon$. Since $\epsilon > 0$ was arbitrarily chosen, this gives the result.

0.2in

[43]{} \[1\][\#1]{} \[1\][`#1`]{} urlstyle \[1\][doi: \#1]{}

E. Aitnouri, S. Wang, and D. Ziou. On comparison of clustering techniques for histogram pdf estimation. *Pattern Recognition and Image Analysis*, 100 (2):0 206–217, 2000.

D. Arthur and S. Vassilvitskii. $k$-means++: The advantages of careful seeding. In *[ACM-SIAM]{} Symposium on Discrete Algorithms (SODA)*, pages 1027–1035, 2007.

A. Azzalini and N. Torelli. Clustering via nonparametric density estimation. *Statistics and Computing*, 170 (1):0 71–80, 2007.

M. Belkin, P. Niyogi, and V. Sindhwani. Manifold regularization: [A]{} geometric framework for learning from labelled and unlabeled examples. *Journal of Machine Learning Research*, 7:0 2399–2434, 2006.

S. Ben-David, T. Lu, D. P[á]{}l, and M. Sot[á]{}kov[á]{}. Learning low-density separators. In D. van Dyk and M. Welling, editors, *International Conference on Artificial Intelligence and Statistics (AISTATS)*, [JMLR]{} Workshop and Conference Proceedings, pages 25–32, 2009.

C. Berge. *Topological Spaces*. Macmillan, New York, 1963.

D. Boley. Principal direction divisive partitioning. *Data Mining and Knowledge Discovery*, 20 (4):0 325–344, 1998.

J. F. Bonnans and A. Shapiro. *Perturbation Analysis of Optimization Problems*. Springer Series in Operations Research. Springer, 2000.

J. V. Burke, A. S. Lewis, and M. L. Overton. Approximating subdifferentials by random sampling of gradients. *Mathematics of Operations Research*, 270 (3):0 567–584, 2002.

J. V. Burke, A. S. Lewis, and M. L. Overton. A robust gradient sampling algorithm for nonsmooth, nonconvex optimization. *[SIAM]{} Journal on Optimization*, 150 (3):0 751–779, 2006.

V. Castelli and T. M. Cover. On the exponential value of labeled samples. *Pattern Recognition Letters*, 160 (1):0 105–111, 1995.

V. Castelli and T. M. Cover. The relative value of labeled and unlabeled samples in pattern recognition with an unknown mixing parameter. *[IEEE]{} Transactions on Information Theory*, 420 (6):0 2102–2117, 1996.

O. Chapelle and A. Zien. Semi-supervised classification by low density separation. In R. G. Cowell and Z. Ghahramani, editors, *International Conference on Artificial Intelligence and Statistics ([AISTATS]{})*, pages 57–64, 2005.

A. Cuevas, M. Febrero, and R. Fraiman. Estimating the number of clusters. *Canadian Journal of Statistics*, 280 (2):0 367–382, 2000.

A. Cuevas, M. Febrero, and R. Fraiman. Cluster analysis: a further approach based on density estimation. *Computational Statistics and Data Analysis*, 360 (4):0 441–459, 2001.

F. E. Curtis and X. Que. An adaptive gradient sampling algorithm for nonsmooth optimization. *Optimization Methods and Software*, 280 (6):0 1302–1324, 2013.

C. Ding and T. Li. Adaptive dimension reduction using discriminant analysis and k-means clustering. In *International Conference on Machine Learning (ICML)*, pages 521–528, 2007.

C. Fowlkes, S. Belongie, F. Chung, and J. Malik. Spectral grouping using the Nystr[ö]{}m method. *[IEEE]{} Transactions on Pattern Analysis and Machine Intelligence*, 260 (2):0 214–225, 2004.

P. Fr[ä]{}nti and O. Virmajoki. Iterative shrinking method for clustering problems. *Pattern Recognition*, 390 (5):0 761–775, 2006.

J. A. Hartigan. *Clustering Algorithms*. Wiley Series in Probability and Mathematical Statistics. Wiley, New York, 1975.

M. Ji, T. Yang, B. Lin, R. Jin, and J. Han. A simple algorithm for semi-supervised learning with improved generalization error bound. In J. Langford and J. Pineau, editors, *International Conference on Machine Learning (ICML)*, pages 1223–1230, 2012.

T. Joachims. Transductive inference for text classification using support vector machines. In *International Conference on Machine Learning (ICML)*, volume 99, pages 200–209, 1999.

A. Lewis and M. Overton. Nonsmooth optimization via quasi-[N]{}ewton methods. *Mathematical Programming*, 141:0 135–163, 2013.

M. Lichman. machine learning repository, 2013. URL <http://archive.ics.uci.edu/ml>.

B. McWilliams, D. Balduzzi, and J. M. Buhmann. Correlated random features for fast semi-supervised learning. In *Advances in Neural Information Processing Systems (NIPS)*, pages 440–448, 2013.

G. Menardi and A. Azzalini. An advancement in clustering via nonparametric density estimation. *Statistics and Computing*, 240 (5):0 753–767, 2014.

V. I. Morariu, B. V. Srinivasan, V. C. Raykar, R. Duraiswami, and L. S. Davis. Automatic online tuning for fast Gaussian summation. In *Advances in Neural Information Processing Systems (NIPS)*, pages 1113–1120, 2008.

A. Ng, M. Jordan, and Y. Weiss. On spectral clustering: analysis and an algorithm. In T. Dietterich, S. Becker, and Z. Ghahramani, editors, *Advances in Neural Information Processing Systems (NIPS)*, volume 14, pages 849–856, 2002.

E. Polak. On the mathematical foundations of nondifferentiable optimization in engineering design. *[SIAM]{} Review*, 290 (1):0 21–89, 1987.

P. Rigollet. Generalization error bounds in semi-supervised classification under the cluster assumption. *Journal of Machine Learning Research*, 8:0 1369–1392, 2007.

A. Rinaldo and L. Wasserman. Generalized density clustering. *The Annals of Statistics*, 380 (5):0 2678–2722, 2010.

A. Rosenberg and J. Hirschberg. -measure: A conditional entropy-based external cluster evaluation measure. In *Empirical Methods in Natural Language Processing and Computational Natural Language Learning*, volume 7, pages 410–420, 2007.

B. W. Silverman. *Density estimation for statistics and data analysis*, volume 26. press, 1986.

A. Singh, R. D. Nowak, and X. Zhu. Unlabeled data: Now it helps, now it doesn’t. In D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou, editors, *Advances in Neural Information Processing Systems (NIPS)*, pages 1513–1520, 2009.

W. Stuetzle and R. Nugent. A generalized single linkage method for estimating the cluster tree of a density. *Journal of Computational and Graphical Statistics*, 190 (2):0 397–418, 2010.

S. K. Tasoulis, D. K. Tasoulis, and V. P. Plagianakos. Enhancing principal direction divisive clustering. *Pattern Recognition*, 430 (10):0 3391–3411, 2010.

S. Tong and D. Koller. Restricted [Bayes]{} optimal classifiers. In *National Conference on Artificial Intelligence ([AAAI]{})*, pages 658–664, 2000.

V. Vapnik and A. Sterin. On structural risk minimization or overall risk in a problem of pattern recognition. *Automation and Remote Control*, 100 (3):0 1495–1503, 1977.

G. Walther. Granulometric smoothing. *The Annals of Statistics*, 250 (6):0 2273–2299, 1997.

P. Wolfe. On the convergence of gradient methods under constraint. *[IBM]{} Journal of Research and Development*, pages 407–411, 1972.

L. Xu, J. Neufeld, B. Larson, and D. Schuurmans. Maximum margin clustering. In L. K. Saul, Y. Weiss, and L. Bottou, editors, *Advances in Neural Information Processing Systems (NIPS)*, volume 17, pages 1537–1544, 2004.

L. Zelnik-Manor and P. Perona. Self-tuning spectral clustering. In *Advances in Neural Information Processing Systems (NIPS)*, pages 1601–1608, 2004.

K. Zhang, I. W. Tsang, and J. T. Kwok. Maximum margin clustering made practical. *[IEEE]{} Transactions on Neural Networks*, 200 (4):0 583–596, 2009.
