Kernel clustering: density biases and solutions

TL;DR

This paper analyzes density biases in kernel clustering methods, proves their existence, and proposes solutions like density equalization to improve clustering accuracy across various applications.

Contribution

It formally demonstrates the presence of density biases in kernel clustering and introduces density equalization techniques as principled solutions.

Findings

Density biases are theoretically proven in kernel K-means.

Density equalization via adaptive weights or kernels mitigates biases.

Experiments confirm the effectiveness of proposed solutions.

Abstract

Kernel methods are popular in clustering due to their generality and discriminating power. However, we show that many kernel clustering criteria have density biases theoretically explaining some practically significant artifacts empirically observed in the past. For example, we provide conditions and formally prove the density mode isolation bias in kernel K-means for a common class of kernels. We call it Breiman's bias due to its similarity to the histogram mode isolation previously discovered by Breiman in decision tree learning with Gini impurity. We also extend our analysis to other popular kernel clustering methods, e.g. average/normalized cut or dominant sets, where density biases can take different forms. For example, splitting isolated points by cut-based criteria is essentially the sparsest subset bias, which is the opposite of the density mode bias. Our findings suggest that a principled solution for density biases in kernel clustering should directly address data inhomogeneity. We show that density equalization can be implicitly achieved using either locally adaptive weights or locally adaptive kernels. Moreover, density equalization makes many popular kernel clustering objectives equivalent. Our synthetic and real data experiments illustrate density biases and proposed solutions. We anticipate that theoretical understanding of kernel clustering limitations and their principled solutions will be important for a broad spectrum of data analysis applications across the disciplines.

Figures

Tables

Table 1. TABLE I : Interactive segmentation errors. AA stands for the average association, NC stands for the normalized cut. Errors are averaged over the GrabCut dataset [ 38 ] , see samples in Figure  8 . ∗ We use [ 37 , 39 ] for a combination of Kernel K-means objective ( 8 ) with Markov Random Field (MRF) regularization terms. The relative weight of the MRF terms is chosen to minimize the average error on the dataset. † Without the MRF term, [ 37 ] and [ 39 ] correspond to the standard kernel K-means  [ 7 , 9 ] .

regularization	average error, %
(boundary smoothness)	Gaussian AA	Gaussian NC	$KNN$ AA	$KNN$ NC
none†	20.4	17.6	12.2	12.4
Euclidean length∗	15.1	16.0	10.2	11.0
contrast-sensitive∗	9.7	13.8	7.1	7.8

Equations

$$k(x,y)\propto \exp \left(-\frac{\Vert x-y{\Vert }^2}{2{\sigma }^2}\right).$$

$$\left(\text{parbox}38.74988pt\text{centering}\mathbf{k}-\mathbf{meanscriterion}\text{@add@centering}\right)\sum _k\sum _{p\in S^k}\Vert f_p-m_k{\Vert }^2.$$

$$F(S,m)=\sum _k\sum _{p\in S^k}\Vert {\varphi }_p-m_k{\Vert }^2$$

$${\hat{m}}_k=\frac{{\sum }_{q\in S^k}{\varphi }_q}{|S^k|}$$

$$F(S)=\sum _k\sum _{p\in S^k}{\Vert {\varphi }_p-\frac{{\sum }_{q\in S^k}{\varphi }_q}{|S^k|}\Vert }^2.$$

$$F(S)\;\;\mbox{$\;\overset{c}{=}\;$}\;\;-\sum_{k}\frac{\sum_{pq\in S^{k}}\langle\phi(f_{p}),\phi(f_{q})\rangle}{|S^{k}|}$$

$$k(x,y):=⟨\varphi (x),\varphi (y)⟩.$$

$$\left(\parbox{38.74988pt}{\centering\bf\small kernel \\ k-means criterion\@add@centering}\right)\quad\quad\quad F(S)\;\;\mbox{$\;\overset{c}{=}\;$}\;\;-\sum_{k}\frac{\sum_{pq\in S^{k}}k(f_{p},f_{q})}{|S^{k}|}.\quad$$

$$-\sum _k\frac{{\sum }_{pq\in S^k}{A}_{pq}}{|S^k|}.$$

$$-\sum _k\frac{{\sum }_{pq\in S^k}{A}_{pq}}{{\sum }_{p\in S^k}{d}_p}.$$

$$\mbox{${\cal P}_{\Sigma}$}(x|S^{k})\;\;:=\;\;\frac{\sum_{q\in S^{k}}k(x,f_{q})}{|S^{k}|},$$

$$k(x,y)\;\;=\;\;|\Sigma|^{-\frac{1}{2}}\;\mbox{$\psi$}\left(\Sigma^{-\frac{1}{2}}(x-y)\right)$$

$$\mbox{$\psi$}(t)\;\propto\;\exp\left(-\frac{\|t\|^{2}}{2\;}\right),$$

$$\sqrt{{\Sigma }_{ii}}=\frac{r_i}{\sqrt[N+4]{n}},{\Sigma }_{ij}=0\text{ for }i\ne j$$

$$F(S)\;\mbox{$\;\overset{c}{=}\;$}\;-\sum_{k}\sum_{p\in S^{k}}\mbox{${\cal P}_{\Sigma}$}(f_{p}|S^{k}).$$

$$-\sum _k\sum _{p\in S^k}\log P(f_p|{\theta }_k).$$

$$\begin{gathered} \left(\text{parbox}38.74988pt\text{centering}\mathbf{entropy} \\ \mathbf{criterion}\text{@add@centering}\right)\sum _k|S^k|\cdot H(S^k) \end{gathered}$$

$$H(S^k):=-\int P(x|{\theta }_k)\log P(x|{\theta }_k)\mathrm{d}x.$$

$$\begin{gathered} \left(\text{parbox}55.97205pt\text{centering}\mathbf{discrete} \\ \mathbf{Ginicriterion}\text{@add@centering}\right)\sum _k|S^k|\cdot G(S^k) \end{gathered}$$

$$G(S^{k})\;:=\;1-\sum_{l\in\cal L}\mbox{${\cal P}$}(l\,|S^{k})^{2}$$

$$\begin{gathered} \left(\text{parbox}55.97205pt\text{centering}\mathbf{continuous} \\ \mathbf{Ginicriterion}\text{@add@centering}\right)\sum _k{w}_k\cdot G(s,k), \end{gathered}$$

$$G(s,k)\;:=\;1-\int\mbox{$\rho^{s}_{k}$}(x)^{2}\,\mathrm{d}x.$$

$$\sum _{p\in C}g(f_p)\approx |C|\int g(x){\rho }_C(x)\mathrm{d}x$$

$$F(S)\;\mbox{$\;\overset{c}{\approx}\;$}\;-\sum_{k}|S^{k}|\int\mbox{${\cal P}_{\Sigma}$}(x|S^{k})\,\mbox{$\rho^{s}_{k}$}(x)\,\,\mathrm{d}x.$$

$$\mbox{${\cal P}_{\Sigma}$}(\cdot\,|S^{k})\;\;\approx\;\;\mbox{$\rho^{s}_{k}$}(\cdot).$$

$$F(S)\;\;\mbox{$\;\overset{c}{\approx}\;$}\;\;-\sum_{k}|S^{k}|\cdot\int\mbox{$\rho^{s}_{k}$}(x)^{2}\,\mathrm{d}x\;\;\mbox{$\;\overset{c}{\equiv}\;$}\;\;\sum_{k}|S^{k}|\cdot G(s,k).$$

$$s_1=\arg \underset{x}{\max }\rho (x)\text{and}{s}_2={\mathcal{R}}^N\setminus s_1.$$

$$\mathbf{E}z:=\int z(x)\rho (x)\mathrm{d}x.$$

$$L(s)\;\;:=\;\;w\int\mbox{$\rho^{s}_{1}$}(x)^{2}\,\mathrm{d}x\;+\;(1-w)\int\mbox{$\rho^{s}_{2}$}(x)^{2}\,\mathrm{d}x$$

$$w:=w_1={\int }_{s(x)=1}\rho (x)\mathrm{d}x=\mathbf{E}[s(x)=1]$$

Full text

1 Introduction

In machine learning, kernel clustering is a well established data analysis technique [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] that can identify non-linearly separable structures, see Figure 1(a-b). Section 1.1 reviews the kernel K-means and related clustering objectives, some of which have theoretically explained biases, see Section 1.2. In particular, Section 1.2.2 describes the discrete Gini clustering criterion standard in decision tree learning where Breiman [11] proved a bias to histogram mode isolation.

Empirically, it is well known that kernel K-means or average association (see Section 1.1.1) has a bias to so-called “tight” clusters for small bandwidths [3]. Figure 1(c) demonstrates this bias on a non-uniform modification of a typical toy example for kernel K-means with common Gaussian kernel

[TABLE]

This paper shows in Section 2 that under certain conditions kernel K-means approximates the continuous generalization of the Gini criterion where we formally prove a mode isolation bias similar to the discrete case analyzed by Breiman. Thus, we refer to the “tight” clusters in kernel K-means as Breiman’s bias.

We propose a density equalization principle directly addressing the cause of Breiman’s bias. First, Section 3 discusses modification of the density with adaptive point weights. Then, Section 4 shows that a general class of locally adaptive geodesic kernels [10] implicitly transforms data and modifies its density. We derive “density laws” relating adaptive weights and kernels to density transformations. They allow to implement density equalization resolving Breiman’s bias, see Figure 1(d). One popular heuristic [12] approximates a special case of our Riemannian kernels.

Besides mode isolation, kernel clustering may have the opposite density bias, e.g. sparse subsets in Normalized Cut [3], see Figure 9(a). Section 5 presents “normalization” as implicit density inversion establishing a formal relation between sparse subsets and Breiman’s bias. Equalization addresses any density biases. Interestingly, density equalization makes many standard kernel clustering criteria conceptually equivalent, see Section 6.

1.1 Kernel K-means

A popular data clustering technique, kernel K-means [1] is a generalization of the basic K-means method. Assuming $\Omega$ denotes a finite set of points and $f_{p}\in\mathcal{R}^{N}$ is a feature (vector) for point $p$, the basic K-means minimizes the sum of squared errors within clusters, that is, distances from points $f_{p}$ in each cluster $S_{k}\subset\Omega$ to the cluster means $m_{k}$

[TABLE]

Instead of clustering data points $\{f_{p}\;|\;p\in\Omega\}\subset\mathcal{R}^{N}$ in their original space, kernel K-means uses mapping $\phi\;\;:\;\;\mathcal{R}^{N}\rightarrow\mathcal{H}$ embedding input data $f_{p}\in\mathcal{R}^{N}$ as points $\phi_{p}\equiv\phi(f_{p})$ in a higher-dimensional Hilbert space $\mathcal{H}$. Kernel K-means minimizes the sum of squared errors in the embedding space corresponding to the following (mixed) objective function

[TABLE]

where $S=(S^{1},S^{2},\dots,S^{K})$ is a partitioning (clustering) of $\Omega$ into $K$ clusters, $m=(m_{1},m_{2},\dots m_{K})$ is a set of parameters for the clusters, and $\|.\|$ denotes the Hilbertian norm111 Our later examples use finite-dimensional embeddings $\phi$ where $\mathcal{H}=\mathcal{R}^{M}$ is an Euclidean space ($M\gg N$) and $\|.\|$ is the Euclidean norm.. Kernel K-means finds clusters separated by hyperplanes in $\mathcal{H}$. In general, these hyperplanes correspond to non-linear surfaces in the original input space $\mathcal{R}^{N}$. In contrast to (3), standard K-means objective (2) is able to identify only linearly separable clusters in $\mathcal{R}^{N}$.

Optimizing $F$ with respect to the parameters yields closed-form solutions corresponding to the cluster means in the embedding space:

[TABLE]

where $|.|$ denotes the cardinality (number of points) in a cluster. Plugging optimal means (4) into objective (3) yields a high-order function, which depends solely on the partition variable $S$:

[TABLE]

Expanding the Euclidean distances in (5), one can obtain an equivalent pairwise clustering criterion expressed solely in terms of inner products $\langle\phi(f_{p}),\phi(f_{q})\rangle$ in the embedding space $\mathcal{H}$:

[TABLE]

where $\;\overset{c}{=}\;$ means equality up to an additive constant. The inner product is often replaced with kernel $k$, a symmetric function:

[TABLE]

Then, kernel K-means objective (5) can be presented as

[TABLE]

Formulation (8) enables optimization in high-dimensional space $\mathcal{H}$ that only uses kernel computation and does not require computing the embedding $\phi(x)$. Given a kernel function, one can use the kernel K-means without knowing the corresponding embedding. However, not any symmetric function corresponds to the inner product in some space. Mercer’s theorem [2] states that any positive semidefinite (p.s.d.) kernel function $k(x,y)$ can be expressed as an inner product in a higher-dimensional space. While p.s.d. is a common assumption for kernels, pairwise clustering objective (8) is often extended beyond p.s.d. affinities. There are many other extension of kernel K-means criterion (8). Despite the connection to density modes made in our paper, kernel clustering has only a weak relation to mean-shift [13], *e.g. *see [14].

1.1.1 Related graph clustering criteria

Positive semidefinite kernel $k(f_{p},f_{q})$ in (8) can be replaced by an arbitrary pairwise similarity or affinity matrix $A=[A_{pq}]$. This yields the average association criterion, which is known in the context of graph clustering [3, 15, 7]:

[TABLE]

The standard kernel K-means algorithm [7, 9] is not guaranteed to decrease (9) for improper (non p.s.d.) kernel $k:=A$. However, [15] showed that dropping p.s.d. assumption is not essential: for arbitrary association $A$ there is a p.s.d. kernel $k$ such that objective (8) is equivalent to (9) up to a constant.

In [3] authors experimentally observed that the average association (9) or kernel K-means (8) objectives have a bias to separate small dense group of data points from the rest, *e.g. *see Figure 2.

Besides average association, there are other pairwise graph clustering criteria related to kernel K-means. Normalized cut is a common objective in the context of spectral clustering [3, 16]. It optimizes the following objective

[TABLE]

where $d_{p}=\sum_{q\in\Omega}A_{pq}$. Note that for $d_{p}=1$ equation (10) reduces to (9). It is known that Normalized cut objective is equivalent to a weighted version of kernel K-means criterion [17, 7].

1.1.2 Probabilistic interpretation via kernel densities

Besides kernel clustering, kernels are also commonly used for probability density estimation. This section relates these two independent problems. Standard multivariate kernel density estimate or Parzen density estimate for the distribution of data points within cluster $S^{k}$ can be expressed as follows [18]:

[TABLE]

with kernel $k$ having the form:

[TABLE]

where $\psi$ is a symmetric multivariate density and $\Sigma$ is a symmetric positive definite bandwidth matrix controlling the density estimator’s smoothness. One standard example is the Gaussian (normal) kernel (1) corresponding to

[TABLE]

which is commonly used both in kernel density estimation [18] and kernel clustering [6, 3].

The choice of bandwidth $\Sigma$ is crucial for accurate density estimation, while the choice of $\psi$ plays only a minor role [19]. There are numerous works regarding kernel selection for accurate density estimation using either fixed [20, 19, 21] or variable bandwidth [22]. For example, Scott’s rule of thumb is

[TABLE]

where $n$ is the number of points, and $r_{i}^{2}$ is the variance of the $i$-th feature that could be interpreted as the range or scale of the data. Scott’s rule gives optimal mean integrated squared error for normal data distribution, but in practice it works well in more general settings. In all cases the optimal bandwidth for sufficiently large datasets is a small fraction of the data range [23, 18]. For shortness, we use adjective $r$-small to describe bandwidths providing accurate density estimation.

If kernel $k$ has form (12) up to a positive multiplicative constant then kernel K-means objective (8) can be expressed in terms of kernel densities (11) for points in each cluster [6]:

[TABLE]

1.2 Other clustering criteria and their known biases

One of the goals of this paper is a theoretical explanation for the bias of kernel K-means with small bandwidths toward tight dense clusters, which we call Breiman’s bias, see Figs 1-2. This bias was observed in the past only empirically. As discussed in Section 4.1, large bandwidth reduces kernel K-means to basic K-means where bias to equal cardinality clusters is known [24]. This section reviews other standard clustering objectives, entropy and Gini criteria, that have biases already well-understood theoretically. In Section 2 we establish a connection between Gini clustering and kernel K-means in case of $r$-small kernels. This connection allows theoretical analysis of Breiman’s bias in kernel K-means.

1.2.1 Probabilistic K-means and entropy criterion

Besides non-parametric kernel K-means clustering there are well-known parametric extensions of basic K-means (2) based on probability models. Probabilistic K-means [24] or model based clustering [25] use some given likelihood functions $P(f_{p}|\theta_{k})$ instead of distances $\|f_{p}-\theta_{k}\|^{2}$ in (2) as in clustering objective

[TABLE]

Note that objective (16) reduces to basic K-means (2) for Gaussian probability model $P(.|\theta_{k})$ with mean $\theta_{k}$ and a fixed scalar covariance matrix.

In probabilistic K-means (16) models can differ from Gaussians depending on a priori assumptions about the data in each cluster, *e.g. *gamma, Gibbs, or other distributions can be used. For more complex data, each cluster can be described by highly-descriptive parametric models such as Gaussian mixtures (GMM). Instead of kernel density estimates in kernel K-means (15), probabilistic K-means (16) uses parametric distribution models. Another difference is the absence of the $\log$ in (15) compared to (16).

The analysis in [24] shows that in case of highly descriptive model $P$, *e.g. *GMM or histograms, (16) can be approximated by the standard entropy criterion for clustering:

[TABLE]

where $H(S^{k})$ is the entropy of the distribution of the data in $S^{k}$:

[TABLE]

The discrete version of the entropy criterion is widely used for learning binary decision trees in classification [11, 18, 26]. It is known that the entropy criterion above is biased toward equal size clusters [11, 24, 27].

1.2.2 Discrete Gini impurity and criterion

Both Gini and entropy clustering criteria are widely used in the context of decision trees [18, 26]. These criteria are used to decide the best split at a given node of a binary classification tree [28]. The Gini criterion can be written for clustering $\{S^{k}\}$ as

[TABLE]

where $G(S^{k})$ is the Gini impurity for the points in $S^{k}$. Assuming discrete feature space $\cal L$ instead of $\mathcal{R}^{N}$, the Gini impurity is

[TABLE]

where \mbox{{\cal P}}(\cdot\,|S^{k}) is the empirical probability (histogram) of discrete-valued features $f_{p}\in\cal L$ in cluster $S^{k}$.

Similarly to the entropy, Gini impurity $G(S^{k})$ can be viewed as a measure of sparsity or “peakedness” of the distribution for points in $S^{k}$. Note that (18) has a form similar to the entropy criterion in (17), except that entropy $H$ is replaced by the Gini impurity. Breiman [11] analyzed the theoretical properties of the discrete Gini criterion (18) when \mbox{{\cal P}}(\cdot\,|S^{k}) are discrete histograms. He proved

Theorem 1** (Breiman).**

For $K=2$ the minimum of the Gini criterion (18) for discrete Gini impurity (19) is achieved by assigning all data points with the highest-probability feature value in $\mathcal{L}$ to one cluster and the remaining data points to the other cluster, as in example for $\mathcal{L}=\{1,\dots,14\}$ on the left. ∎

2 Breiman’s bias (numerical features)

In this section we show that the kernel K-means objective reduces to a novel continuous Gini criterion under some general conditions on the kernel function, see Section 2.1. We formally prove in Section 2.2 that the optimum of the continuous Gini criterion isolates the data density mode. That is, we show that the discussed earlier biases observed in the context of clustering [3] and decision tree learning [11] are the same phenomena. Section 2.3 establishes connection to maximum cliques [29] and dominant sets [8].

For further analysis we reformulate the problem of clustering a discrete set of points $\{f_{p}\,|\,p\in\Omega\}\subset\mathcal{R}^{N}$, see Section 1.1, as a continuous domain clustering problem. Let $P$ be a probability measure over domain $\mathcal{R}^{N}$ and $\rho$ be the corresponding continuous probability density function such that the discrete points $f_{p}$ could be treated as samples from this distribution. The clustering of the continuous domain will be described by an assignment function $s:\mathcal{R}^{N}\to\{1,2,\dots,K\}$. Density $\rho$ implies conditional probability densities \mbox{\rho^{s}_{k}}(x):=\rho(x\,|\,s(x)=k). Feature points $f_{p}$ in cluster $S^{k}$ could be interpreted as a sample from conditional density $\rho^{s}_{k}$.

Then, the continuous clustering problem is to find an assignment function optimizing a clustering criteria. For example, we can analogously to (18) define continuous Gini clustering criterion

[TABLE]

where $w_{k}$ is the probability to draw a point from $k$-th cluster and

[TABLE]

In the next section we show that kernel K-means energy (15) can be approximated by continuous Gini-clustering criterion (20) for $r$-small kernels.

2.1 Kernel K-means and continuous Gini criterion

To establish the connection between kernel clustering and the Gini criterion, let us first recall Monte-Carlo estimation [24], which yields the following expectation-based approximation for a continuous function $g(x)$ and cluster $C\subset\Omega$:

[TABLE]

where $\rho_{C}$ is the “true” continuous density of features in cluster $C$. Using (22) for $C=S^{k}$ and g(x)=\mbox{{\cal P}_{\Sigma}}(x|S^{k}), we can approximate the kernel density formulation in (15) by its expectation

[TABLE]

Note that partition $S=(S^{1},\dots,S^{K})$ is determined by dataset $\Omega$ and assignment function $s$. We also assume

[TABLE]

This is essentially an assumption on kernel bandwidth. That is, we assume that kernel bandwidth gives accurate density estimation. For shortness, we call such bandwidths $r$-small, see Section 1.1.2. Then (23) reduces to approximation

[TABLE]

Additional application of Monte-Carlo estimation $|S^{k}|/|\Omega|\approx w_{k}$ allows replacing set cardinality $|S^{k}|$ by probability $w_{k}$ of drawing a point from $S^{k}$. This results in continuous Gini clustering criterion (20), which approximates (15) or (8) up to an additive and positive multiplicative constants.

Next section proves that the continuous Gini criterion (20) has a similar bias observed by Breiman in the discrete case.

2.2 Breiman’s bias in continuous Gini criterion

This section extends Theorem 1 to continuous Gini criterion (20). Since Section 2.1 has already established a close relation between continuous Gini criterion and kernel K-means for $r$-small bandwidth kernels, then Breiman’s bias also applies to the latter. For simplicity, we focus on $K=2$ as in Breiman’s Theorem 1.

Theorem 2** (Breiman’s bias in continuous case).**

For $K=2$ the continuous Gini clustering criterion (20) achieves its optimal value at the partitioning of $\mathcal{R}^{N}$ into regions

[TABLE]

Proof.

The statement follows from Lemma 2 below. ∎

We denote mathematical expectation of function $z\;:\;\Omega\to{\mathcal{R}}^{1}$

[TABLE]

Minimization of (20) corresponds to maximization of the following objective function

[TABLE]

where the probability to draw a point from cluster $1$ is

[TABLE]

where $[\cdot]$ is the indicator function. Note that mixed joint density

[TABLE]

allows to write conditional density $\rho^{s}_{1}$ in (26) as

[TABLE]

Equations (26) and (27) give

[TABLE]

Introducing notation

[TABLE]

allows to further rewrite objective function $L(s)$ as

[TABLE]

Without loss of generality assume that $\frac{{\mathbf{E}}(1-I)\rho}{1-{\mathbf{E}}I}\leq\frac{{\mathbf{E}}I\rho}{{\mathbf{E}}I}$ (the opposite case would yield a similar result). We now need following

Lemma 1**.**

Let $a,b,c,d$ be some positive numbers, then

[TABLE]

Proof.

Use reduction to a common denominator. ∎

Lemma 1 implies inequality

[TABLE]

which is needed to prove the Lemma below.

Lemma 2**.**

Assume that function $s_{\varepsilon}$ is

[TABLE]

Then

[TABLE]

Proof.

Due to monotonicity of expectation we have

[TABLE]

Then (30) and (33) imply

[TABLE]

That is, the right part of (32) is an upper bound for $L(s)$.

Let $I_{\varepsilon}\equiv[s_{\varepsilon}(x)=1]$. It is easy to check that

[TABLE]

Definition (31) also implies

[TABLE]

This result and (33) conclude that

[TABLE]

Finally, the limits in (35) and (37) imply

[TABLE]

This equality and bound (34) prove (32). ∎

This result states that the optimal assignment function separates the mode of the density function from the rest of the data. The proof considers case $K=2$ for continuous Gini criterion approximating kernel K-means for $r$-small kernels. The multi-cluster version for $K>2$ also has Breiman’s bias. Indeed, it is easy to show that any two clusters in the optimal solution shall give optimum of objective (20). Then, these two clusters are also subject to Breiman’s bias. See a multi-cluster example in Figure 3.

Practical considerations: While Theorem 2 suggests that the isolated density mode should be a single point, in practice Breiman’s bias in kernel k-means isolates a slightly wider cluster around the mode, see Figures 2, 3, 7(a-d), 8. Indeed, Breiman’s bias holds for kernel k-means when the assumptions in Section 2.1 are valid. In practice, shrinking of the clusters invalidates approximations (23) and (24) preventing the collapse of the clusters.

2.3 * Connection to maximal cliques and dominant sets*

Interestingly, there is also a relation between maximum cliques and density modes. Assume [math]-​$1$ kernel $[\|x-y\|\leq\sigma]$ with bandwidth $\sigma$. Then, kernel matrix $A$ is a connectivity matrix corresponding to a $\sigma$-disk graph. Intuitively, the maximum clique on this graph should be inside a disk with the largest number of points in it, which corresponds to the density mode.

Formally, mode isolation bias can be linked to both maximum clique and its weighted-graph generalization, dominant set [8]. It is known that maximum clique [29] and dominant set [8] solve a two-region clustering problem with energy

[TABLE]

corresponding to average association (9) for $K=1$ and $S^{1}\subseteq\Omega$. Under the same assumptions as above, Gini impurity (21) can be used as an approximation reducing objective (39) to

[TABLE]

Using (33) and (37) we can conclude that the optimum of (40) isolates the mode of density function $\rho$. Thus, clustering minimizing (39) for $r$-small bandwidths also has Breiman’s bias. That is, for such bandwidths the concepts of maximum clique and dominant set for graphs correspond to the concept of mode isolation for data densities. Dominant sets for the examples in Figures 1(c), 2(a), and 7(d) would be similar to the shown mode-isolating solutions.

3 Adaptive weights solving Breiman’s bias

We can use a simple modification of average association by introducing weights $w_{p}\geq 0$ for each point “error” within the equivalent kernel K-means objective (3)

[TABLE]

Such weighting is common for K-means [23]. Similarly to Section 1.1 we can expand the Euclidean distances in (41) to obtain an equivalent weighted average association criterion generalizing (9)

[TABLE]

Weights $w_{p}$ have an obvious interpretation based on (41); they change the data by replicating each point $p$ by a number of points in the same location (Figure 4a) in proportion to $w_{p}$. Therefore, this weighted formulation directly modifies the data density as

[TABLE]

where $\rho_{p}$ and $\rho^{\prime}_{p}$ are respectively the densities of the original and the new (replicated) points. The choice of $w_{p}=1/\rho_{p}$ is a simple way for equalizing data density to solve Breiman’s bias. As shown in Figure 4(a), such a choice enables low-density points to be replicated more frequently than high-density ones. This is one of density equalization approaches giving the solution in Figure 1(d).

4 Adaptive kernels solving Breiman’s bias

Breiman’s bias in kernel K-means is specific to $r$-small bandwidths. Thus, it has direct implications for the bandwidth selection problem discussed in this section. Note that kernel bandwidth selection for clustering should not be confused with kernel bandwidth selection for density estimation, an entirely different problem outlined in Section 1.1.2. In fact, $r$-small bandwidths give accurate density estimation, but yield poor clustering due to Breiman’s bias. Larger bandwidths can avoid this bias in clustering. However, Section 4.1 shows that for extremely large bandwidths kernel K-means reduces to standard K-means, which loses ability of non-linear cluster separation and has a different bias to equal cardinality clusters [24, 27].

In practice, avoiding extreme bandwidths is problematic since the notions of small and large strongly depend on data properties that may significantly vary across the domain, e.g. in Figure 1c,d where no fixed bandwidth gives a reasonable separation. This motivates locally adaptive strategies. Interestingly, Section 4.2 shows that any locally adaptive bandwidth strategy implicitly corresponds to some data embedding $\Omega\to\mathcal{R}^{N^{\prime}}$ deforming density of the points. That is, locally adaptive selection of bandwidth is equivalent to selection of density transformation. Local kernel bandwidth and transformed density are related via the density law established in (59). As we already know from Theorem 2, Breiman’s bias is caused by high non-uniformity of the data, which can be addressed by density equalizing transformations. Section 4.3 proposes adaptive kernel strategies based on our density law and motivated by a density equalization principle addressing Breiman’s bias. In fact, a popular locally adaptive kernel in [12] is a special case of our density equalization principle.

4.1 Overview of extreme bandwidth cases

Section 2.1 and Theorem 2 prove that for $r$-small bandwidths the kernel K-means is biased toward “tight” clusters, as illustrated in Figures 1, 2 and 7(d). As bandwidth increases, continuous kernel density (11) no longer approximates the true distribution $\rho^{s}_{k}$ violating (24). Thus, Gini criterion (25) is no longer valid as an approximation for kernel K-means objective (15). In practice, Breiman’s bias disappears gradually as bandwidth gets larger. This is also consistent with experimental comparison of smaller and larger bandwidths in [3].

The other extreme case of bandwidth for kernel K-means comes from its reduction to basic K-means for large kernels. For simplicity, assume Gaussian kernels (1) of large bandwidth $\sigma$ approaching data diameter. Then the kernel can be approximated by its Taylor expansion $\exp\left(-\frac{\|x-y\|^{2}}{2\sigma^{2}}\right)\approx 1-\frac{\|x-y\|^{2}}{2\sigma^{2}}$ and kernel K-means objective (8) for $\sigma\gg\|x-y\|$ becomes222Relation (44) easily follows by substituting $m_{k}\equiv\frac{1}{|S^{k}|}\sum_{p\in S^{k}}f_{p}$. (up to a constant)

[TABLE]

which is equivalent to basic K-means (2) for any fixed $\sigma$.

Figure 5 summarizes kernel K-means biases for different bandwidths. For large bandwidths the kernel K-means loses its ability to find non-linear cluster separation due to reduction to the basic K-means. Moreover, it inherits the bias to equal cardinality clusters, which is well-known for the basic K-means [24, 27]. On the other hand, for small bandwidths kernel K-means has Breiman’s bias proven in Section 2. To avoid the biases in Figure 5, kernel K-means should use a bandwidth neither too small nor too large. This motivates locally adaptive bandwidths.

4.2 Adaptive kernels as density transformation

This section shows that kernel clustering (8) with any locally adaptive bandwidth strategy satisfying some reasonable assumptions is equivalent to fixed bandwidth kernel clustering in a new feature space (Theorem 3) with a deformed point density. The adaptive bandwidths relate to density transformations via density law (59). To derive it, we interpret adaptiveness as non-uniform variation of distances across the feature space. In particular, we use a general concept of geodesic kernel defining adaptiveness via a metric tensor and illustrate it by simple practical examples.

Our analysis of Breiman’s bias in Section 2 applies to general kernels (12) suitable for density estimation. Here we focus on clustering with kernels based on radial basis functions $\psi$ s.t.

[TABLE]

To obtain adaptive kernels, we replace Euclidean metric with Riemannian inside (45). In particular, $\|x-y\|$ is replaced with geodesic distances $d_{g}(x,y)$ between features $x,y\in{\cal R}^{N}$ based on any given metric tensor $g(f)$ for $f\in{\cal R}^{N}$. This allows to define a geodesic or Riemannian kernel at any points $f_{p}$ and $f_{q}$ as in [10]

[TABLE]

where $d_{pq}:=d_{g}(f_{p},f_{q})$ is introduced for shortness.

In practice, the metric tensor can be defined only at the data points $g_{p}:=g(f_{p})$ for $p\in\Omega$. Often, quickly decaying radial basis functions $\psi$ allow Mahalanobis distance approximation inside (46)

[TABLE]

which is normally valid only in a small neighborhood of $f_{p}$. If necessary, one can use more accurate approximations for $d_{g}(f_{p},f_{q})$ based on Dijkstra [33] or Fast Marching method [34].

Example 1** **(**Adaptive non-normalized333 Lack of normalization as in (48)

is critical for density equalization resolving Breiman’s bias, which is our only goal for adaptive kernels. Note that without kernel normalization as in (12) Parzen density formulation of kernel k-means (15) no longer holds invalidating the relation to Gini and Breiman’s bias in Section 2. On the contrary, normalized variable kernels are appropriate for density estimation [22] validating (15). They can also make approximation (24) more accurate strengthening connections to Gini and Breiman’s bias. Gaussian kernel**).

Mahalanobis distances based on (adaptive) bandwidth matrices $\Sigma_{p}$ defined at each point $p$ can be used to define adaptive kernel

[TABLE]

which equals fixed bandwidth Gaussian kernel (1) for $\Sigma_{p}=\sigma^{2}I$. Kernel (48) approximates (46) for exponential function $\psi$ in (13) and tensor $g$ continuously extending matrices $\Sigma^{-1}_{p}$ over the whole feature space so that $g_{p}=\Sigma^{-1}_{p}$ for $p\in\Omega$. Indeed, assuming matrices $\Sigma^{-1}_{p}$ and tensor $g$ change slowly between points within bandwidth neighbourhoods, one can use (47) for all points in

[TABLE]

due to exponential decay outside the bandwidth neighbourhoods.

Example 2** (Zelnik-Manor & Perona kernel [12]).**

This popular kernel is defined as $\kappa_{pq}:=\exp\frac{-\|f_{p}-f_{q}\|^{2}}{2\sigma_{p}\sigma_{q}}$. This kernel’s relation to (46) is less intuitive due to the lack of “local” Riemannian tensor. However, under assumptions similar to those in (49), it can still be seen as an approximation of geodesic kernel (46) for some tensor $g$ such that $g_{p}=\sigma_{p}^{-2}I$ for $p\in\Omega$. They use heuristic $\sigma_{p}=R_{p}^{K}$, which is the distance to the K-th nearest neighbour of $f_{p}$.

Example 3** (KNN kernel).**

*This adaptive kernel is defined as u_{p}(f_{p},f_{q})=[f_{q}\in\mbox{K!N!N}(f_{p})] where \mbox{K!N!N}(f_{p}) is the set of $K$ nearest neighbors of $f_{p}$. This kernel approximates (46) for uniform function $\psi(t)=[t<1]$ and tensor $g$ such that $g_{p}=I/(R^{K}_{p})^{2}$. *

Theorem 3**.**

Clustering (8) with (adaptive) geodesic kernel (46) is equivalent to clustering with fixed bandwidth kernel $k^{\prime}(f^{\prime}_{p},f^{\prime}_{q}):=\psi^{\prime}(\|f^{\prime}_{p}-f^{\prime}_{q}\|)$ in new feature space $\mathcal{R}^{N^{\prime}}$ for some radial basis function $\psi^{\prime}$ using the Euclidean distance and some constant $N^{\prime}$.

Proof.

A powerful general result in [35, 36, 15] states that for any symmetric matrix $(d_{pq})$ with zeros on the diagonal there is a constant $h$ such that squared distances

[TABLE]

form Euclidean matrix $(\widetilde{d}_{pq})$. That is, there exists some Euclidean embedding $\Omega\to{\cal R}^{N^{\prime}}$ where for $\forall p\in\Omega$ there corresponds a point $f^{\prime}_{p}\in{\cal R}^{N^{\prime}}$ such that $\|f^{\prime}_{p}-f^{\prime}_{q}\|=\widetilde{d}_{pq}$, see Figure 6. Therefore,

[TABLE]

for $\psi^{\prime}(t)\!:=\!\psi(\sqrt{t^{2}-h^{2}[t\geq h]})$ and $k_{g}(f_{p},f_{q})\!=\!k^{\prime}(f^{\prime}_{p},f^{\prime}_{q})$. ∎

Theorem 3 proves that adaptive kernels for $\{f_{p}\}\subset{\cal R}^{N}$ can be equivalently replaced by a fixed bandwidth kernel for some implicit embedding444The implicit embedding implied by Euclidean matrix (50) should not be confused with embedding in the Mercer’s theorem for kernel methods. $\{f^{\prime}_{p}\}\subset{\cal R}^{N^{\prime}}$ in a new space. Below we establish a relation between three local properties at point $p\,$: adaptive bandwidth represented by matrix $g_{p}$ and two densities $\rho_{p}$ and $\rho^{\prime}_{p}$ in the original and the new feature spaces. For $\varepsilon>0$ consider an ellipsoid in the original space $\mathcal{R}^{N}$, see Figure 6(a),

[TABLE]

Assuming $\varepsilon$ is small enough so that approximation (47) holds, ellipsoid (52) covers features $\{f_{q}\,|\,q\in\Omega_{p}\}$ for subset of points

[TABLE]

Similarly, consider a ball in the new space $\mathcal{R}^{N^{\prime}}$, see Figure 6(b),

[TABLE]

covering features $\{f^{\prime}_{q}\,|\,q\in\Omega^{\prime}_{p}\}$ for points

[TABLE]

It is easy to see that (50) implies $\Omega_{p}=\Omega^{\prime}_{p}$. Let $\rho_{p}$ and $\rho^{\prime}_{p}$ be the densities555We use the physical rather than probability density. They differ by a factor. of points within $B_{p}$ and $B^{\prime}_{p}$ correspondingly. Assuming $|\cdot|$ denotes volumes or cardinalities of sets, we have

[TABLE]

Omitting a constant factor depending on $\varepsilon$, $h$, $N$ and $N^{\prime}$ we get

[TABLE]

representing the general form of the density law. For the basic isotropic metric tensor such that $g_{p}=I/\sigma_{p}^{2}$ it simplifies to

[TABLE]

Thus, bandwidth $\sigma_{p}$ can be selected adaptively based on any desired transformation of density $\rho^{\prime}_{p}\equiv\tau(\rho_{p})$ using

[TABLE]

where observed density $\rho_{p}$ in the original feature space can be evaluated at any point $p$ using any standard estimators, *e.g. *(11).

4.3 Density equalizing locally adaptive kernels

Bandwidth formula (59) works for any density transform $\tau$. To address Breiman’s bias, one can use density equalizing transforms \tau(\rho)=\mbox{\mathit{const}} or $\tau(\rho)=\frac{1}{\alpha}\log(1+\alpha\rho)$, which even up

the highly dense parts of the feature space as illustrated on the right. Some empirical results using density equalization \tau(\rho)=\mbox{\mathit{const}} for synthetic and real data are shown in Figures 1(d) and 7(e,f).

One way to estimate the density in (59) is $K\!N\!N$ approach [18]

[TABLE]

where $n\equiv|\Omega|$ is the size of the dataset, $R_{p}^{K}$ is the distance to the $K$-th nearest neighbor of $f_{p}$, $V_{K}$ is the volume of a ball of radius $R_{p}^{K}$ centered at $f_{p}$. Then, density law (59) for \tau(\rho)=\mbox{\mathit{const}} gives

[TABLE]

consistent with heuristic bandwidth in [12], see Example 2.

The result in Figure 1(d) uses adaptive Gaussian kernel (48) for $\Sigma_{p}=\sigma_{p}I$ with $\sigma_{p}$ derived in (61). Theorem 3 claims equivalence to a fixed bandwidth kernel in some transformed higher-dimensional space $\mathcal{R}^{N^{\prime}}\!$. Bandwidths (61) are chosen specifically to equalize the data density in this space so that \tau(\rho)=\mbox{\mathit{const}}.

The picture on the right illustrates such density equalization for the data in Figure 1(d). It shows a 3D projection of the transformed data obtained by multi-dimensional scaling [32] for matrix $(\widetilde{d}_{pq})$ in (50). The observed density equalization removes Breiman’s bias from the clustering in Figure 1(d).

Real data experiments for kernels with adaptive bandwidth (61) are reported in Figures 2, 3, 7, 8 and Table I. Figure 7(e) illustrates the empirical density equalization effect for this bandwidth. Such data homogenization removes the conditions leading to Breiman’s bias, see Theorem 2. Also, we observe empirically that $K\!N\!N$ kernel is competitive with adaptive Gaussian kernels, but its sparsity gives efficiency and simplicity of implementation.

5 Normalized Cut and Breiman’s bias

Breiman’s bias for kernel K-means criterion (8), a.k.a. average association (AA) (9), was empirically identified in [3], but our Theorem 2 is its first theoretical explanation. This bias was the main critique against AA in [3]. They also criticize graph cut [40] that “favors cutting small sets of isolated nodes”. These critiques are used to motivate normalized cut (NC) criterion (10) aiming at balanced clustering without “clumping” or “splitting”.

We do not obeserve any evidence of the mode isolation bias in NC. However, Section 5.1 demonstrates that NC still has a bias to isolating sparse subsets. Moreover, using the general density analysis approach introduced in Section 4.2 we also show in Section 5.2 that normalization implicitly corresponds to some density-inverting embedding of the data. Thus, mode isolation (Breiman’s bias) in this implicit embedding corresponds to the sparse subset bias of NC in the original data.

5.1 Sparse subset bias in Normalized Cut

The normalization in NC does not fully remove the bias to small isolated subsets and it is easy to find examples of “splitting” for weakly connected nodes, see Figure 9(a). The motivation argument for the NC objective below Fig.1 in [3] implicitly assumes similarity matrices with zero diagonal, which excludes many common similarities like Gaussian kernel (1). Moreover, their argument is built specifically for an example with a single isolated point, while an isolated pair of points will have a near-zero NC cost even for zero diagonal similarities.

Intuitively, this NC issue can be interpreted as a bias to the “sparsest” subset (Figure 9a), the opposite of AA’s bias to the “densest” subset, *i.e. *Breiman’s bias (Figure 1c). The next subsection discusses the relation between these opposite biases in detail. In any case, both of these density inhomogeneity problems in NC and AA are directly addressed by our density equalization principle embodied in adaptive weights $w_{p}\propto 1/\rho_{p}$ in Section 3 or in the locally adaptive kernels derived in Section 4.3. Indeed, the result in Figure 1(d) can be replicated with NC using such adaptive kernel. Interestingly, [12] observed another data non-homogeneity problem in NC different from the sparse subset bias in Figure 9(a), but suggested a similar adaptive kernel as a heuristic solving it.

5.2 Normalization as density inversion

The bias to sparse clusters in NC with small bandwidths (Figure 9a) seems the opposite of mode isolation in AA (Figure 1c). Here we show that this observation is not a coincidence since NC can be reduced to AA after some density-inverting data transformation. While it is known [17, 7] that NC is equivalent to weighted kernel K-means (i.e. weighted AA) with some modified affinity, this section relates such kernel modification to an implicit density-inverting embedding where mode isolation (Breiman’s bias) corresponds to sparse clusters in the original data.

First, consider standard weighted AA objective for any given affinity/kernel matrix $\hat{A}_{pq}=k(f_{p},f_{q})$ as in (42)

[TABLE]

Clearly, weights based on node degrees $w=d$ and “normalized” affinities $\hat{A}_{pq}=\frac{A_{pq}}{d_{p}d_{q}}$ turn this into NC objective (10). Thus, average association (9) becomes NC (10) after two modifications:

• •

replacing $A_{pq}$ by normalized affinities $\hat{A}_{pq}=\frac{A_{pq}}{d_{p}d_{q}}$ and

• •

introducing point weights $w_{p}=d_{p}$.

Both of these modifications of AA can be presented as implicit data transformations modifying denisty. In particular, we show that the first one “inverses” density turning sparser regions into denser ones, see Figure 10(a). The second data modification is generally discussed as a density transform in (43). We show that node degree weights $w_{p}=d_{p}$ do not remove the “density inversion”.

For simplicity, assume standard Gaussian kernel (1) based on Euclidean distances $d_{pq}=\|f_{p}-f_{q}\|$ in $\mathcal{R}^{N}$

[TABLE]

To convert AA into NC we first need an affinity “normalization”

[TABLE]

equivalently formulated as a modification of distances

[TABLE]

Using a general approach in the proof of Theorem 3, there exists some Euclidean embedding $\bar{f}_{p}\in\mathcal{R}^{\bar{N}}$ and constant $h\geq 0$ such that

[TABLE]

Thus, modified affinities $\hat{A}_{pq}$ in (62) correspond to the Gaussian kernel for the new embedding $\{\bar{f}_{p}\}$ in $\mathcal{R}^{\bar{N}}$

[TABLE]

Assuming $d_{q}\approx d_{p}$ for features $f_{q}$ near $f_{p}$, equations (63) and (64) imply the following relation for such neighbors of $f_{p}$

[TABLE]

Then, similarly to the arguments in (56), a small ball of radius $\varepsilon$ centered at $f_{p}$ in $\mathcal{R}^{N}$ and a ball of radius $\sqrt{\varepsilon^{2}+h^{2}+4\sigma^{2}\log(d_{p})}$ at $\bar{f}_{p}$ in $\mathcal{R}^{\bar{N}}$ contain the same number of points. Thus, similarly to (57) we get a relation between densities at points $f_{p}$ and $\bar{f}_{p}$

[TABLE]

This implicit density transformation is shown in Figure 10(a). Sub-linearity in dense regions addresses mode isolation (Breiman’s bias). However, sparser regions become relatively dense and kernel-modified AA may split them. Indeed, the result in Figure 9(a) can be obtained by AA with normalized affinity $\frac{A_{pq}}{d_{p}d_{q}}$.

The second required modification of AA introduces point weights $w_{p}=d_{p}$. It has an obvious equivalent formulation via data points replication discussed in Section 3, see Figure 4(a). Following (43), we obtain its implicit density modification effect $\rho^{\prime}_{p}=d_{p}\bar{\rho}_{p}$. Combining this with density transformation (65) implied by affinity normalization $\frac{A_{pq}}{d_{p}d_{q}}$, we obtain the following density transformation effect corresponding to NC, see Figure 10(b),

[TABLE]

The density inversion in sparse regions relates NC’s result in Figure 9(a) to Breiman’s bias for embedding $\{\bar{f}_{p}\}$ in $\mathcal{R}^{\bar{N}}$.

Figure 10 shows representative plots for density transformations (65), (66) using the following node degree approximation based on Parzen approach (11) for Gaussian affinity (kernel) $A$

[TABLE]

Empirical relation between $d_{p}$ and $\rho_{p}$ is illistrated below: some

overestimation occurs for sparcer regions and underestimation happens for denser regions. The node degree for Gaussian kernels has to be at least $1$ (for an isolated node) and at most $N$ (for a dense graph).

6 Discussion (kernel clustering equivalence)

Density equalization with adaptive weights in Section 3 or adaptive kernels in Section 4 are useful for either AA or NC due to their density biases (mode isolation or sparse subset). Interestingly, kernel clustering criteria discussed in [3] such as normalized cut (NC), average cut (AC), average association (AA) or kernel K-means are practically equivalent for such adaptive methods. This can be seen both empirically (Table I) and conceptually. Note, weights $w_{p}\propto 1/\rho_{p}$ in Section 3 produce modified data with near constant node degrees $d^{\prime}_{p}\propto\rho^{\prime}_{p}\propto 1$, see (67) and (43). Alternatively, KNN kernel (Example 3) with density equalizing bandwidth (61) also produce nearly constant node degrees $d_{p}\approx K$ where $K$ is the neighborhood size. Therefore, both cases give

[TABLE]

which correspond to NC (10), AA (9), and AC criteria. As discussed in [3], the last objective also has very close relations with standard partitioning concepts in spectral graph theory: isoperimetric or Cheeger number, Cheeger set, ratio cut.

This equivalence argument applies to the corresponding clustering objectives and is independent of specific optimization algorithms developed for them. Interestingly, the relation between (9) and basic K-means objective (3) suggests that standard Lloyd’s algorithm can be used as a basic iterative approach for approximate optimization of all clustering criteria in (68). In practice, however, kernel K-means algorithm corresponding to the exact high-dimensional embedding $\{\phi_{p}\}$ in (3) is more sensitive to local minima compared to iterative K-means over approximate lower-dimensional embeddings based on PCA [14, Section 3.1]666K-means is also commonly used as a discretization heuristic for spectral relaxation [3] where a similar eigen analysis is motivated by spectral graph theory [41, 42, 43] defferently from PCA dimensionalty reduction in [14]..

7 Conclusions

This paper identifies and proves density biases, *i.e. *isolation of modes or sparsest subsets, in many well-known kernel clustering criteria such as kernel K-means (average association), ratio cut, normalized cut, dominant sets. In particular, we show conditions when such biases happen. Moreover, we propose density equalization as a general principle for resolving such biases. We suggest two types of density equalization techniques using adaptive weights or adaptive kernels. We also show that density equalization unifies many popular kernel clustering objectives by making them equivalent.

Acknowledgements

The authors would like to thank Professor Kaleem Siddiqi (McGill University) for suggesting a potential link between Breiman’s bias and the dominant sets. This work was generously supported by the Discovery and RTI programs of the National Science and Engineering Research Council of Canada (NSERC).