Semi-supervised model-based clustering with controlled clusters leakage

TL;DR

This paper introduces C3L, a semi-supervised Gaussian mixture model that incorporates user-defined leakage control to effectively discover natural subgroups in partially labeled data, enhancing clustering accuracy.

Contribution

The paper presents a novel semi-supervised clustering model with leakage control, along with theoretical analysis and an efficient optimization algorithm.

Findings

C3L effectively identifies meaningful clusters in partially labeled data.

The model's leakage parameter allows flexible control over clustering consistency.

Experimental results demonstrate high-quality clustering performance.

Abstract

In this paper, we focus on finding clusters in partially categorized data sets. We propose a semi-supervised version of Gaussian mixture model, called C3L, which retrieves natural subgroups of given categories. In contrast to other semi-supervised models, C3L is parametrized by user-defined leakage level, which controls maximal inconsistency between initial categorization and resulting clustering. Our method can be implemented as a module in practical expert systems to detect clusters, which combine expert knowledge with true distribution of data. Moreover, it can be used for improving the results of less flexible clustering techniques, such as projection pursuit clustering. The paper presents extensive theoretical analysis of the model and fast algorithm for its efficient optimization. Experimental results show that C3L finds high quality clustering model, which can be applied in discovering meaningful groups in partially classified data.

Figures

Tables

Table 1. Table 1: Regression UCI data sets and median number of clusters returned by C3L .

	Airfoil	Forest	Music+	Stock
# Instances	1502	517	1059	536
# Features	5	12	5	7
# Clusters	5	3	5	7

Table 2. Table 2: UCI data sets used in subgroups detection experiment. Last two rows show which reference groups were used for creating classes Y − subscript 𝑌 Y_{-} and Y + subscript 𝑌 Y_{+} .

	Balance	Segmentation+	User	Wine
# Instances	625	210	258	178
# Features	4	5	5	13
# Clusters	3	7	4	3
$Y_{-}$	{1,2}	{1,3,4,5,7}	{1,2}	{1,2}
$Y_{+}$	{3}	{2,6}	{3,4}	{3}

Table 3. Table 3: Normalized mutual information for UCI datasets.

Method	Balance	Segmentation	User	Wine
C3L0.01	0.50	0.62	0.53	0.50
C3L0.05	0.44	0.60	0.49	0.51
CECH	0.48	0.59	0.47	0.34
CEC	0.03	0.58	0.25	0.46
c-GMM	0.20	0.54	0.56	0.47
GMMH	0.49	0.56	0.66	0.40
GMM	0.07	0.58	0.36	0.45

Table 4. Table 4: Correlation between the accuracy of decision boundary and clustering results.

Method	Balance	Segmentation	User	Wine
C3L0.01	0.95	0.64	0.98	0.45
C3L0.05	0.90	0.34	0.94	0.36
CECH	0.99	0.68	0.96	0.76
GMMH	0.99	0.98	0.88	0.80

Equations

$$H=\{x\in {\mathbb{R}}^N:h^{T}x=a\},$$

$${\mathrm{class}}_H(x)=\mathrm{sign}(h^{T}x-a).$$

$$X_{\pm }=\{x\in X:{\mathrm{class}}_H(x)=\pm 1\}.$$

$$H_{\pm }=\{x\in {\mathbb{R}}^N:{\mathrm{class}}_H(x)=\pm 1\},$$

$${\alpha }_i=\min \{{\int }_{H_{-}}{g}_i(x)dx,{\int }_{H_{+}}{g}_i(x)dx\}.$$

$$h^{\times }(X\Vert g)=-\frac{1}{|X|}\sum _{x\in X}\ln (g(x)).$$

$$\begin{array}[]{l}h^{\times}(X\|\mathcal{N}(m,\Sigma))=\\ \tfrac{N}{2}\mathrm{ln}(2\pi)+\tfrac{1}{2}\|\hat{m}_{X}-m\|_{\Sigma}+\tfrac{1}{2}\mathrm{tr}(\Sigma^{-1}\SS{X})+\tfrac{1}{2}\mathrm{ln}\mathrm{det}(\Sigma),\end{array}$$

$$\begin{array}[]{l}\hat{m}_{X}:=\frac{1}{|X|}\sum\limits_{x\in X}x,\\ \SS{X}:=\frac{1}{|X|}\sum\limits_{x\in X}(x-\hat{m}_{X})(x-\hat{m}_{X})^{T},\end{array}$$

$$\Vert x{\Vert }_{\Sigma }:=x^T{\Sigma }^{-1}x$$

$$\sum _i{p}_i(-\ln p_i)+p_i{h}^{\times }(X_i\Vert g_i),$$

$$\text{either }{\int }_{H_{-}}g(x)dx\ge 1-\alpha \text{ or }{\int }_{H_{+}}g(x)dx\ge 1-\alpha .$$

$$\begin{array}[]{l}E^{\alpha}_{H}(X\|\mathcal{G}):=\\ \inf\{h^{\times}(X\|g):\,g\in\mathcal{G}\text{ which satisfies constraint }(H,\alpha)\},\end{array}$$

$$\begin{array}[]{l}E^{\alpha}_{H}(X_{1},\ldots,X_{k}\|\mathcal{G})=\\ \sum_{i=1}^{k}p_{i}(-\mathrm{ln}p_{i}+E^{\alpha}_{H}(X_{i}\|\mathcal{G}))\text{, where }p_{i}=\frac{|X_{i}|}{|X|}.\end{array}$$

$$\mathrm{class}(x)=\arg \underset{j=\pm }{\max }{f}_j(x).$$

$$\mathrm{class}(x)=\mathrm{sign}f(x).$$

$$X\ni x\to (f(x),x)\in \mathbb{R}\times {\mathbb{R}}^N.$$

$${\mathcal{G}}_{1,N}:=\{g(x)=g_1(x_1)\cdot g_N(x_{2:N+1})\text{ : }{g}_1\in {\mathcal{G}}_1,g_N\in {\mathcal{G}}_N\},$$

$$\begin{array}[]{l}\int_{[0,+\infty)\times\mathbb{R}^{N-1}}g(x)dx=\int_{0}^{+\infty}g_{1}(x_{1})dx_{1},\\ \int_{(-\infty,0]\times\mathbb{R}^{N-1}}g(x)dx=\int_{-\infty}^{0}g_{1}(x_{1})dx_{1}.\end{array}$$

$$h^{\times }(X\Vert g)=h^{\times }(X[1]\Vert g_1)+h^{\times }(X[2:N]\Vert g_{N-1}),$$

$$\begin{array}[]{l}h^{\times}(X\|g)\\ =\sum\limits_{x\in X}-\mathrm{ln}g(x)=\sum\limits_{x\in X}-\mathrm{ln}(g_{1}(x_{1})\cdot g_{N-1}(x_{2:N}))\\[4.30554pt] =\sum\limits_{x\in X}\big{(}-\mathrm{ln}(g_{1}(x_{1}))-\mathrm{ln}(g_{N-1}(x_{2:N}))\big{)}\\[4.30554pt] =\sum\limits_{x_{1}\in X{[1]}}-\mathrm{ln}(g_{1}(x_{1}))+\sum\limits_{x_{N-1}\in X{[2:N]}}-\mathrm{ln}(g_{N-1}(x_{N-1}))\\[4.30554pt] =h^{\times}(X{[1]}\|g_{1})+h^{\times}(X{[2:N]}\|g_{N-1}).\end{array}$$

$$E_H^{\alpha }(X\Vert {\mathcal{G}}_{1,N-1})=E_{\{0\}}^{\alpha }(X[1]\Vert {\mathcal{G}}_1)+h^{\times }(X[2:N]\Vert {\mathcal{G}}_{N-1}),$$

$$p^{\alpha }:={\Phi }_{\mathcal{N}(0,1)}^{-1}(1-\alpha ),$$

$$\alpha ={\Phi }_{\mathcal{N}(m,\sigma )}(m-p^{\alpha }\sigma ),$$

$$|m|\ge p^{\alpha }\sigma ,$$

$$\begin{array}[]{l}h^{\times}(X\|\mathcal{N}(m,\sigma))=\frac{1}{2}\left(\frac{\hat{\sigma}_{X}^{2}+(m-\hat{m}_{X})^{2}}{\sigma^{2}}+\mathrm{ln}(\sigma^{2})+\mathrm{ln}(2\pi)\right).\end{array}$$

$$\begin{array}[]{ll}m_{X}^{\alpha}:=&\frac{-(p^{\alpha})^{2}\hat{m}_{X}+\mathrm{sign}(\hat{m}_{X})p^{\alpha}\sqrt{((p^{\alpha})^{2}+4)\hat{m}_{X}^{2}+4\hat{\sigma}_{X}^{2}}}{2},\\ \sigma_{X}^{\alpha}:=&\frac{|\hat{m}_{X}^{\alpha}|}{p^{\alpha}}.\end{array}$$

$$h(m,\sigma )=h^{\times }(X\Vert \mathcal{N}(m,\sigma ))\text{ under the condition }|m|\ge p^{\alpha }\sigma .$$

$$h(m)=\frac{1}{2}\big{(}\frac{\hat{\sigma}_{X}^{2}+(m-\hat{m}_{X})^{2}}{m^{2}}(p^{\alpha})^{2}+\mathrm{ln}\frac{m^{2}}{(p^{\alpha})^{2}}+\mathrm{ln}(2\pi)\big{)}.$$

$$h^{\prime }(m)=-\frac{{\hat{\sigma }}_X^2+{\hat{m}}_X^2}{m^3}(p^{\alpha }{)}^2+\frac{{\hat{m}}_X}{m^2}(p^{\alpha }{)}^2+\frac{1}{m}.$$

$$m=\frac{-(p^{\alpha }{)}^2{\hat{m}}_X+\mathrm{sign}({\hat{m}}_X)p^{\alpha }\sqrt{((p^{\alpha }{)}^2+4){\hat{m}}_X^2+4{\hat{\sigma }}_X^2}}{2},$$

Full text

1 Introduction

Model-based clustering aims at finding a mixture of probability models, which optimally estimates true probability distribution on data space. Contrary to other clustering techniques, it does not only recover meaningful groups, but also gives a rule (probability model) for generating elements from clusters. Therefore, it is commonly used in various areas of machine learning and data analysis [40, 28, 34].

Although clustering is an unsupervised technique, one can introduce additional information to guide the algorithm what is the expected structure of clusters. Semi-supervised learning methods usually use partial labeling [16] or pairwise constraints [20] to transfer expert knowledge into clustering process, while consensus and alternative clustering gather information from several partitions of data into one general view [22, 9]. In this paper, we assume that we have the knowledge about division of data set into two categories and focus on the following problem: How to find the best model of clusters that preserves a fixed amount of information about existing categories? In other words, we focus on finding interesting clusters, which are very likely to belong to one category.

To explain a basic motivation behind our model, let us consider an expert system used for automatic text translation. It is a common practice to construct several translation models, each designed for one cluster retrieved from a data set [1]. Alternatively, since texts are often categorized into specific domains, e.g. sport, politics, etc., then each translator can be fitted to one of these categories. To consider together both options, we could implement a separate module responsible for finding clusters, which (a) are described by compact models (e.g. Gaussians) and (b) are related with predefined topics. Observe that optimization of these two conflicting goals simultaneously is non-trivial. We cannot cluster elements from each category individually, because this strategy does not lead to optimal solution for the entire data set (in terms of likelihood). Moreover, existing categorization might be inaccurate as well as the interesting groups can cross the boundary between predefined domains. Therefore, a better approach is to incorporate the constraint to the clustering process and always work with the entire data set.

Our method can also be applied to strictly unsupervised situations, where no initial categorization is given. Let us recall that one way to analyze clusters in complex data spaces relies on finding projections onto one dimensional subspaces, where groups can be easily identified. Projection pursuit focuses on choosing such a direction, which optimizes selected statistical index such as kurtosis [26] or skewness [17]. Since one dimensional views generate linear decision boundaries in original data space, it is not possible to find flexible cluster structures. However, we can input such a linear boundary to our model in order to improve existing clusters. Our method directly uses the information from initial splitting, but can extend linear decision surfaces to nonlinear ones generated by probabilistic mixture models.

Following the above motivation, we propose a semi-supervised clustering with controlled clusters leakage model (C3L), which integrates a distribution of data with a fixed division of the space into two categories. C3L focuses on finding a type of Gaussian mixture model (GMM) [21], which maximizes the likelihood function and preserves the information contained in the initial splitting with a predefined probability (leakage level). Intuitively, we allow for the flow of clusters densities over decision surface, but with a full control of total probability assigned to the opposite category, which is defined as the leakage level $\alpha\in(0,1)$ (see Figure 1). This general idea is formulated as a constrained optimization problem (Section 3).

The advantages of C3L can be summarized as follows:

It has a closed form solution in a special case of cross-entropy clustering (a type of GMM) [35]. 2. 2.

It can be efficiently implemented and optimized by a modified on-line Hartigan algorithm (Section 4). 3. 3.

The user can directly parametrize C3L by a maximal inconsistency level between initial categorization and final clustering model (leakage level). 4. 4.

The selection of the leakage level $\alpha$ allows to move from a strictly unsupervised GMM for $\alpha=0.5$, where decision boundary has no effect on clustering, to the limiting case of $\alpha\to 0$, where every group is fully condensed in one category (Section 5).

Experimental studies confirm that the proposed approach builds a high quality model under a given constraint in terms of inner clustering measures, such as Bayesian Information Criterion (Section 6.1). It can be successfully used to discover meaningful groups in partially classified data (Section 6.2) as well as to improve existing clusters obtained by applying projection techniques (Section 6.3). We present a real-life case study, in which the use of C3L allows to detect subgroups of chemical space given their division into active and inactive classes (Section 6.4).

2 Related work

Semi-supervised clustering incorporates the knowledge about class labels to partitioning process [5]. This information can be presented as partial labeling, which gives a division of a small portion of data into categories, or as pairwise constraints, which indicate whether two data points originate from the same (must-links) or distinct classes (cannot-links). Although pairwise constraints provide less amount of information than partial labeling, it is easier to assess whether two instances come from the same group than assign them to particular classes.

Clustering with pairwise constraints was introduced by Wagstaff et al. [37], who created a variant of k-means, which focuses on preserving all constraints. Shental et al. [29] constructed a version of Gaussian mixture model, which gathers data points into equivalence classes (called chunklets) using must-link relation and then applied EM algorithm on such generalized data set of chunklets. This approach was later modified to multi-modal clustering models [32]. The aforementioned methods work well with noiseless side information, but deteriorate the results when some constraints are mislabeled. To overcome this problem, the authors of [6, 19] applied hidden Markov random fields (HMRF) to construct more sophisticated dependencies between linked points. However, the use of HMRF leads to complex solutions, which are difficult to optimize. In recent years, Asafi and Cohen-Or. [4] suggested reducing distances between data points with a must-link constraint and adding a dimension for each cannot-link constraint. After updating all other distances to, e.g., satisfy the triangle inequality, the thus obtained pairwise distance matrix can be used for unsupervised learning. Wang and Davidson[38] proposed a version of spectral clustering, which relies on solving a generalized eigenvalue problem.

Partial labeling is used in clustering to define sample data points from particular classes. Liu and Fu [16] added additional attributes to feature vectors and proposed modified k-means algorithm. There is also a semi-supervised version of fuzzy c-means [25, 24], where the authors supplied the cost function with a regularization term that penalizes fuzzy partitions that are inconsistent with the side information. GMMs can be adapted to make use of class labels by combining the classical unsupervised GMM with a supervised one [2, 41].

Since assigning data points to classes or labeling pairwise constraints requires extensive domain knowledge, then many clustering methods were adapted to use additional information about data, which does not require human intervention. One example is consensus clustering, which considers gathering information coming from different domains [22]. On the other hand, complementary (alternative) clustering aims at finding groups which provide a perspective on the data that expands on what can be inferred from previous partitions [9].

C3L is a version of Gaussian mixture model, which uses side information given by class labels or more generally by a decision boundary between classes. In contrast to classical methods applying partial labeling, it focuses on finding subgroups of original classes. This goal is similar to information bottleneck method [36, 7]. Roughly speaking, this approach tries to construct compact clusters (compressed representation), which contain high amount of information about existing classes (auxiliary variable). While information-theoretic approaches use mutual information (or conditional entropy) to preserve the consistency with an initial categorization, C3L explicitly defines maximal probability of inconsistency (leakage level). The leakage level can be understood as Bayes error in classification or significance level in hypothesis testing and restricts every cluster model to be assigned to one of two initial classes with a predefined probability. Subgroups could also be detected by using cannot-link constraints, clustering with pairwise constraints does not allow to input maximal level of error. Moreover, computational complexity of applying cannot-link constraints to GMM is usually high, while C3L works in a comparable time to classical unsupervised mixture models.

C3L can be naturally combined with projection pursuit approach, which focuses on selecting low dimensional projections of data for finding clusters (or other meaningful characteristics of data). Such a projection can be determined by optimizing selected statistical coefficients e.g. kurtosis [26, 11] or skewness [17, 18]. Since every one dimensional view induces linear decision boundary in the original space, this technique may not be sufficient to detect complex data patters. C3L allows to take such a rough linear splitting of data and correct simple decision boundaries to nonlinear ones.

3 Theoretical model

In this section, we introduce our model and discuss its possible extensions and applications.

Let a data space $\mathbb{R}^{N}$ be divided by a codimension one hyperplane111$(N-1)$-dimensional hyperplane $H$ given by:

[TABLE]

for fixed $h\in\mathbb{R}^{N}$ and $a\in\mathbb{R}$. The hyperplane $H$ induces hard classification rule: the class label of each point $x\in\mathbb{R}^{N}$ is determined by

[TABLE]

This splits a dataset $X\subset\mathbb{R}^{N}$ into two groups $X_{+},X_{-}$ given by

[TABLE]

Alternatively, we can consider the initial classification of entire space $\mathbb{R}^{N}$ as

[TABLE]

and put $X_{\pm}=X\cap H_{\pm}$. Note that the uncertainty of class label is usually higher for elements localized closer to the barrier than for those with larger distance from $H$.

In a model-based clustering we focus on estimating a density of data space with a use of mixture of $k$ densities, $g=\sum_{i=1}^{k}p_{i}g_{i}$, where every $g_{i}$ belongs to a given parametric family of densities (usually Gaussian) and $p_{i}$ are prior probabilities [21]. This goal can be practically realized by maximizing the likelihood function. Let us define the inconsistency between a cluster density $g_{i}$ and the initial classification:

[TABLE]

The above formula gives the amount of probability that is spread to opposite class and it is related with Bayes error of Gaussian model assuming that $H$ predicts the class membership correctly.

Our question is: how to find a clustering model that optimizes a likelihood function and provides high consistency with initial classification? If we knew that $H$ gives a perfect classification rule, then we could try to keep every model $g_{i}$ maximally consistent with $H$, i.e. perform a separate clustering of every category. However, this is usually not the case and this strategy does not guarantee to obtain optimal solution (in terms of likelihood) for the entire data set. Moreover, some interesting groups can cross the decision boundary. Therefore, we should allow for the flow of corresponding densities over a decision surface, but with the full control of the total probability assigned to the opposite class. In our approach, we formulate a constrained optimization problem, where we aim at finding such a mixture model $g$ that maximizes the quality of density estimation and preserves a fixed inconsistency level $\alpha$, i.e. every component $g_{i}$ has to satisfy $\alpha_{i}\leq\alpha$.

One could probably try to realize the above goal by a classical GMM approach, however, at very high numerical and theoretical cost. It would be impossible to get a closed form solution and complex non-linear optimization would be needed. Therefore, in this paper we have decided to use cross-entropy clustering (CEC) [35, 33, 30], which similarly to GMM divides data with respect to Gaussian distributions. Contrary to GMM, in CEC the clusters do not “cooperate” one with another to build the global cost function222Instead of optimizing a density $p_{i}g_{1}+\ldots+p_{k}g_{k}$, CEC finds optimal subdensity $\max\{p_{i}g_{1},\ldots,p_{k}g_{k}\}$, i.e. every point $x$ is linked with exactly one model $g_{i}$. and consequently it is enough to calculate the cost function for each cluster individually.

CEC is based on Minimum Description Length Principle (MDLP) [27] and focuses on minimizing the generalized cross-entropy function, by selecting optimal Gaussian probability distribution for each cluster333The minimization of cross-entropy is equivalent to the maximization of likelihood function.. Given a single cluster $X$ and corresponding density $g$, the empirical cross-entropy function equals:

[TABLE]

In the case of Gaussian densities, $g=\mathcal{N}(m,\Sigma)$, the above formula can be reduced to its closed form:

[TABLE]

where

[TABLE]

denote the sample mean and covariance of $X$ and

[TABLE]

is the Mahalanobis norm. Given $k$ clusters $X_{1},\ldots,X_{k}$ the overall cross-entropy function equals

[TABLE]

where $g_{i}$ is a Gaussian density with parameters given by (4) for $X_{i}$ and $p_{i}=\frac{|X_{i}|}{|X|}$.

The term $p_{i}(-\mathrm{ln}p_{i})$ adds a cost for maintaining a cluster. In consequence, the method tends to keep the model simple and allows for the reduction of redundant clusters. Therefore, CEC cost function (5) combines the model accuracy with its complexity.

With this, we are ready to define our C3L model. First, we give a definition of a linear constraint, which restricts every cluster model to one category with a fixed probability.

Definition 1**.**

Let $H$ be a codimension one hyperplane on $X\subset\mathbb{R}^{N}$ and let $\alpha>0$. We say that a density $g$ satisfies a linear constraint $(H,\alpha)$, if

[TABLE]

The number $\alpha$ will be referred as the leakage level.

The above definition of linear constraint is analogical to linear separability with power $\alpha$ given in [26, Section 2] in the context of projection pursuit.

If $\alpha\geq\frac{1}{2}$ then an arbitrary density satisfies one of the conditions given by (6). Therefore, a strict constraint is given by $\alpha<\frac{1}{2}$. Observe that the above constraint is reminiscent of a typical approach used in hypothesis testing, where we accept a given hypothesis if it lies within a predefined percentage of density. In our case we consider only those density cluster models, which lie with a probability $(1-\alpha)$ on one side of decision boundary.

We now introduce a linear constraint to CEC framework. First, we define the criterion function for a single cluster:

Definition 2**.**

(One cluster C3L cost function) Let $X\subset\mathbb{R}^{N}$ be divided by a codimension one hyperplane $H$. Given $\alpha>0$ and a family $\mathcal{G}$ of Gaussian densities, C3L cost function for $X$ is defined by

[TABLE]

where $h^{\times}(X\|g)$ is given by (3).

We always assume that a covariance matrix of $X$ is nonsingular, i.e., $\mathrm{det}(\SS{X})\neq 0$. This prevents from the situation when $X$ lies in the subspace of $\mathbb{R}^{N}$, which might lead to degenerate solutions.

The overall C3L cost is defined as follows:

Definition 3**.**

(Overall C3L cost function for clustering) Let $X\subset\mathbb{R}^{N}$ be divided by a codimension one hyperplane $H$. Given a splitting of $X$ into clusters $X_{1},\ldots,X_{k}$, a family $\mathcal{G}$ of Gaussian densities and $\alpha>0$, the total C3L clustering cost equals

[TABLE]

The optimal clustering is the one, which minimizes the above cost function (the optimization problem will be the subject of the next section). We emphasize that C3L accepts arbitrary Gaussian densities as mixture components, which satisfy the linear constraint. In particular, each Gaussian can have distinct covariance matrix.

Let us observe that introduced model can be applied in the case of any decision boundary (not only linear hyperplane). If $f_{+}$ and $f_{-}$ are two decision functions that quantify the chance of assigning data points to positive and negative classes, then the label of an instance $x\in\mathbb{R}^{N}$ is chosen as

[TABLE]

Clearly, for two class problem we can define one discriminant $f=f_{+}-f_{-}$. In consequence, the formula (9) can be simplified to:

[TABLE]

Next, we extend the input space $\mathbb{R}^{N}$ to $\mathbb{R}\times\mathbb{R}^{N}$ and embed our data set $X$ into this space by:

[TABLE]

Observe that a hyperplane $H=\{0\}\times\mathbb{R}^{N}$ gives the same classification rule in $\mathbb{R}^{N+1}$ to the formula (10) in $\mathbb{R}^{N}$.

Let $\mathcal{G}_{N}$ denotes the set of all Gaussian densities on $\mathbb{R}^{N}$ and let $\mathcal{G}_{1,N}$ be the set of Gaussian densities, that can be factorized into two independent components, defined by:

[TABLE]

where $x_{k:l}=(x_{k},\ldots,x_{l})$, for $x=(x_{1},\ldots,x_{N+1})$. If we consider a mixture model $g=\sum_{i}p_{i}g_{i}$, where $g_{i}=g^{i}_{1}\cdot g^{i}_{N}\in\mathcal{G}_{1,N}$, then the first component $g^{i}_{1}$ will describe a distribution of discrimination function $f$ while $g^{i}_{N}$ will describe a density in the original space $X$. Therefore, the use of $\mathcal{G}_{1,N}$ allows to model a distribution of data and discriminant function individually. We will use the family $\mathcal{G}_{1,N}$ in the next section.

4 Optimization

Without loss of generality, we assume that a decision boundary of $X\subset\mathbb{R}^{N}$ is given by a hyperplane444Observe that given an arbitrary hyperplane one can always shift the original data and change the basis of $\mathbb{R}^{N}$ in an orthonormal way to obtain this situation. $H=\{0\}\times\mathbb{R}^{N-1}$. From now on, our attention is restricted to the class of Gaussian densities $\mathcal{G}_{1,N-1}$ defined by (11), i.e. we assume that every component $g$ is of the form $g=g_{1}\cdot g_{N-1}$, where $g_{1}$ is $1$-dimensional and $g_{N-1}$ is $(N-1)$-dimensional Gaussian density. This model suits perfectly to the case of arbitrary decision boundary described at the end of section 3, where we model the data distribution and the values of decision support function separately.

We are going to show that in the case of $\mathcal{G}_{1,N-1}$ we can compute one cluster cost function analytically. Let us first observe that the selection of 1-dimensional density $g_{1}\in\mathcal{G}_{1}$ and $(N-1)$-dimensional density $g_{N-1}\in\mathcal{G}_{N-1}$ can be done separately in C3L clustering. First, we verify that the linear constraint is independent of $g_{N-1}$, i.e.,

[TABLE]

Then, observe that the cross-entropy between $X$ and $g$ can be calculated as:

[TABLE]

where $X{[k:l]}$ denotes a data set $X$ restricted to the attributes from $k$ to $l$. This follows from

[TABLE]

Since the product densities can be selected individually, the optimization subject to the constraint is performed in one dimension.

Corollary 1**.**

Let $(H,\alpha)$ be a linear constraint defined on dataset $X\subset\mathbb{R}^{N}$, where $H=\{0\}\times\mathbb{R}^{N-1}$. Then the one cluster C3L cost function of $X$ is given by:

[TABLE]

where $E^{\alpha}_{\{0\}}((X{[1]}\|\mathcal{G}_{1})$ is one cluster C3L cost function (7) calculated in one dimensional situation.

To complete the formula given in Corollary 1, the C3L cost function in one dimensional case has to be calculated. To facilitate the calculation we first give an equivalent form of linear constraint.

Given $\alpha>0$, let us denote by $p^{\alpha}$ the corresponding quantile:

[TABLE]

where $\Phi_{\mathcal{N}(m,\sigma)}(\cdot)$ denotes a cumulative distribution function of $\mathcal{N}(m,\sigma)$. Making use of elementary calculations we get that:

[TABLE]

for any $m\in\mathbb{R}$ and $\sigma>0$. Then, one dimensional density $\mathcal{N}(m,\sigma)$ satisfies the constraint $(\{0\},\alpha)$, iff

[TABLE]

In other words, the distance between the mean $m$ and the barrier is at least $p^{\alpha}\sigma$.

To calculate the optimal one dimensional cost function, we must observe that the cross-entropy (3) between a distribution of dataset $X\subset\mathbb{R}$ with mean $\hat{m}_{X}$ and standard deviation $\hat{\sigma}_{X}$ and a Gaussian density $\mathcal{N}(m,\sigma)$ equals:

[TABLE]

The optimal parameters of $\mathcal{N}(m,\sigma)$ are obtained by the minimization of the above function under the restriction (13):

Theorem 1**.**

Let $X\subset\mathbb{R}$ be a dataset with the mean $\hat{m}_{X}\neq 0$ and the standard deviation $\hat{\sigma}_{X}>0$. We assume that $(\{0\},\alpha)$ denotes the linear constraint on $X$, for $\alpha>0$, and $p^{\alpha}=\Phi^{-1}_{\mathcal{N}(0,1)}(1-\alpha)$.

If $|\hat{m}_{X}|\geq p^{\alpha}\hat{\sigma}_{X}$, then put $m_{X}^{\alpha}:=\hat{m}_{X}$, $\sigma_{X}^{\alpha}:=\hat{\sigma}_{X}$, otherwise

[TABLE]

Then, the normal density $\mathcal{N}(m_{X}^{\alpha},\sigma_{X}^{\alpha})$ minimizes the value of $h^{\times}(X\|\mathcal{N}(m,\sigma))$, given by (14), under the restriction $|m|\geq p^{\alpha}\sigma_{X}$ (C3L cost function $E^{\alpha}_{\{0\}}(X\|\mathcal{G})$).

Proof.

Our aim is to find the minimum of the function

[TABLE]

It is obvious that the above function has the derivative zero only at its global minimum which is given by a pair $(m_{X},\sigma_{X})$. Consequently, if $m=\hat{m}_{X},\sigma=\hat{\sigma}_{X}$ satisfies the constraint $|m|\geq p^{\alpha}\sigma$, then we have found the minimum.

In the opposite case, we only need to verify what happens on the boundary of the constraints (since we do not have any local minimum inside $|m|>p^{\alpha}\sigma$), that is when $p^{\alpha}\sigma=|m|$. Then, by (14), the function (16) simplifies to

[TABLE]

Then

[TABLE]

Finally, the solution of $h^{\prime}(m)=0$ which minimizes the value of $h(m)$, is given by

[TABLE]

which completes the proof. ∎

The above analysis shows how to calculate the best model of clusters for a given partition. However, finding an optimal partition is NP-hard problem, where heuristic iterative algorithms are commonly used. One can apply a slight modification of Hartigan approach to optimize C3L cost function (see A for details). Similar algorithm is used in optimization of CEC and k-means methods.

5 Theoretical analysis

In this section we present a theoretical analysis of C3L model in its simplified form. We start with investigating the convergence of cluster parameters with respect to the leakage level $\alpha$. Then, we show that under certain assumptions a decision boundary determined by C3L model with two clusters converges to the initial barrier, when $\alpha$ approaches to [math].

In order to accommodate the constraint one-dimensional density cluster model modifies its mean and standard deviation according to Theorem 1. The relation between $m_{X}^{\alpha}$ and $\sigma_{X}^{\alpha}$ (given by (15)) is inversely proportional, i.e., the increase of $m_{X}^{\alpha}$ results in the decrease of $\sigma_{X}^{\alpha}$ and vice versa (see Figure 2). However, the most important fact is that $|m_{X}^{\alpha}|$ does not grow infinitely, but converges to a finite number dependent on a data set. To prove it formally, let us first consider a one dimensional case.

Lemma 1**.**

We assume that $X\subset\mathbb{R}$ is a data set with a mean $\hat{m}_{X}\neq 0$ and standard deviation $\hat{\sigma}_{X}>0$. Let $g^{\alpha}=\mathcal{N}(m_{X}^{\alpha},\sigma_{X}^{\alpha})$ denote a density minimizing one cluster C3L cost function under the linear constraint $(\{0\},\alpha)$, i.e, $m_{X}^{\alpha}$ and $\sigma_{X}^{\alpha}$ are given by Theorem 1.

Then:

[TABLE]

Proof.

The proof is included in B. ∎

If we combine the above result with the fact that the mean and the covariance of $(N-1)$ dimensional density $g_{N-1}$, for a model $g=g_{1}\cdot g_{N-1}$, do not depend on the linear constraint, but are the maximum likelihood estimators of data, we get the following corollary:

Corollary 2**.**

We assume that $X\subset\mathbb{R}^{N}$ is a data set, with a mean $\hat{m}_{X}$ and a covariance matrix $\SS{X}$, where $\hat{m}_{X{[1]}}\neq 0$ and $\hat{\sigma}_{X{[1]}}^{2}=\SS{X{[1]}}$. Let $g^{\alpha}=g_{1}^{\alpha}\cdot g_{N-1}\in\mathcal{G}_{1,N-1}$ denote a density minimizing one cluster C3L cost function under the linear constraint $(\{0\}\times\mathbb{R}^{N-1},\alpha)$, i.e., $g^{\alpha}_{1}=\mathcal{N}(m_{X{[1]}}^{\alpha},\sigma_{X{[1]}}^{\alpha})$ is given by Theorem 1 and $g_{N-1}=\mathcal{N}(m_{X{[2:N]}}^{\alpha},\Sigma_{X{[2:N]}}^{\alpha})$.

Then:

[TABLE]

The Figure 3 shows the influence of the change of the parameter $\alpha$ on the form of resulting density function.

We now discuss the relations between an initial decision boundary defined by a hyperplane $H$ and a splitting determined by C3L model. For a simplicity, we consider the case of only two clusters. We assume that the hyperplane $H=\{0\}\times\mathbb{R}^{N-1}$ defines the classification rule (1):

[TABLE]

which splits a data set $X\subset\mathbb{R}^{N}$ into two subsets $X_{\pm}$ (2).

Let the leakage level $\alpha$ be fixed. In a simplified form of C3L model with only two clusters, we assume that density models $p_{\pm}g^{\alpha}_{\pm}$ are chosen so that to maximize one cluster cost functions of $X_{\pm}$, respectively (not the overall cost of clustering). In the other words, a class density is selected ignoring the influence of objects which belong to the opposite class. Then any incoming object $x\in X$ can be classified to one of two clusters by calculating

[TABLE]

We will show that a decision boundary determined by such model converges to $H$, as the leakage level $\alpha$ approaches to 0, i.e.,

[TABLE]

which is illustrated in Figure 4:

Theorem 2**.**

We assume that $X\subset\mathbb{R}^{N}$ is a data set and $H=\{0\}\times\mathbb{R}^{N-1}$ defines a hyperplane in $\mathbb{R}^{N}$ dividing $X$ into two classes $X_{-},X_{+}$, where $X_{\pm}\neq\emptyset$. Let $g^{\alpha}_{\pm}\in\mathcal{G}_{1,N-1}$ denote two densities minimizing one cluster C3L functions of $X_{\pm}$ under the constraint $(H,\alpha)$, respectively.

Then, there exists a constant $C>0$ such that

[TABLE]

for every $x\in\mathbb{R}^{N}$ satisfying $dist(x,H)\leq C$.

The above convergence holds only for instances, which are localized within the margin of size $C$ around the barrier $H$. This is a natural situation occurring in every Gaussian discrimination. To prove the above theorem we first consider one dimensional situation, where an exact value of constant $C$ will be given.

Lemma 2**.**

We assume that $X_{\pm}\subset\mathbb{R}_{\pm}$ are two non empty sets with means $\hat{m}_{\pm}$ and standard deviations $\hat{\sigma}_{\pm}$. Let $g^{\alpha}_{\pm}\in\mathcal{G}_{1}$ denote two densities minimizing one cluster C3L cost functions of $X_{\pm}$ under the linear constraint $(\{0\},\alpha)$, respectively.

Then,

[TABLE]

for every $0<x\leq 2\left(\hat{m}_{+}+\frac{\hat{\sigma}_{+}^{2}}{\hat{m}_{+}}\right)$.

Proof.

The proof is included in C. ∎

In one dimensional case the constant $C$ from Theorem 2 equals $C=2\left(\hat{m}_{+}+\frac{\hat{\sigma}_{+}^{2}}{\hat{m}_{+}}\right)$. Below we complete the proof of our main result:

Proof.

(of Theorem 2) Let $x\in\mathbb{R}^{N}$ be such that $0<x_{1}<2\left(\hat{m}_{+}+\frac{\hat{\sigma}_{+}^{2}}{\hat{m}_{+}}\right)$, where $\hat{m}_{+},\hat{\sigma}_{+}$ are the mean and standard deviation of $X_{+}{[1]}$. In other words, we assume that $x$ lies at the right side of a decision boundary. We assume that the optimal C3L densities $g^{\alpha}_{\pm}\in\mathcal{G}_{1,N-1}$ for $X_{\pm}$ equal

[TABLE]

with the priors $p_{\pm}$. We will show that

[TABLE]

for sufficiently small $\alpha>0$.

Since

[TABLE]

the formula (17) can be rewritten as

[TABLE]

Making use of Lemma 2, we have

[TABLE]

Because

[TABLE]

the LHS of (18) can be arbitrary large, when $\alpha\to 0$, which completes the proof. ∎

6 Experiments

We tested our method on sample examples retrieved from UCI repository [15] and one real data set of chemical compounds [39]. We verified the quality of the clustering model and demonstrated that C3L can be useful in discovering natural subgroups given a partial knowledge about two class division. We also used C3L to extend linear boundary between clusters obtained by projection pursuit technique to nonlinear one. We also show its application on real data set of chemical compounds. We compared its performance with related model-based clustering techniques.

6.1 Quality of the model

In this experiment we consider a scenario, where every instance is assigned to one of two classes based on the value of a fixed decision support function. C3L builds a clustering model which preserves the information of class membership in a sense that every cluster density belongs to one of two classes with a probability greater than $(1-\alpha)$ (6). We want to verify the quality of such model and compare it with the results produced by related model-based clustering techniques, which however do not allow for a direct specification of the leakage level.

To compare the quality of clustering models we applied Bayesian Information Criterion (BIC), which is a standard criterion for model selection [8]. The lower the BIC is the better the model is.

We used four regression UCI examples, which are summarized in Table 1. A dependent (output) variable was treated as a decision support function, which determines a decision boundary $H$. More precisely, if $X\times Y$ is a data set, where $X\subset\mathbb{R}^{N}$ contains explanatory variables and $Y\subset\mathbb{R}$ includes dependent variable, then a decision boundary $H$ is defined by $H=X\times\{\mathrm{median}(Y)\}$, where $\mathrm{median}(Y)$ is the median of attribute $Y$.

The effects of C3L were compared with those obtained by classical GMM method, which ignores the presence of existing decision boundary. To introduce a decision boundary to the model, we also considered the second variant of GMM (which is referred to as GMMH): given a linear hyperplane $H$, which divides a data set $X$ into two regions $X_{-},X_{+}$ (see (2)), GMMH is defined as follows:

• •

GMM is applied to $X_{-}$ and $X_{+}$ separately, which give two models $g_{\pm}=p_{1}^{\pm}g_{1}^{\pm}+\ldots+p_{k}^{\pm}g_{k}^{\pm}$.

• •

These models are combined into a single one by $g=\frac{|X_{-}|}{|X|}g_{-}+\frac{|X_{+}|}{|X|}g_{+}$.

• •

Finally, every point $x\in X$ is assigned to the most probable cluster by calculating $\frac{|X_{\pm}|}{|X|}p_{i}^{\pm}g_{i}^{\pm}(x)$.

Analogical strategies were also applied to CEC. Both variants of GMM use general Gaussian densities to model clusters distributions, while CEC-based methods use the same densities as C3L method, i.e. densities from the family $\mathcal{G}_{1,N-1}$.

We ran each method on $X\times Y$ with a decision boundary $H$. To investigate the influence of the leakage level on the clustering effects of C3L, six leakage levels were considered, $\alpha\in\{0.01,0.05,0.15,0.25,0.35,0.5\}$. Since C3L and CECH internally find the final number of clusters, we ran them with 10 groups, while other methods used the median number of groups returned by C3L calculated over these six leakage levels.

The results presented in the Figure 5 prove that the quality of C3L model improves as the leakage level is increased. This is a natural behavior, because lower values of $\alpha$ indicate higher restrictions on the clusters models. Since other methods do not control the inconsistency with classification, we measured their resulting leakage levels and marked returned BIC values. One can observe that in most cases CECH and GMMH gave worse BIC than C3L method for corresponding leakage levels. It follows from the fact that C3L optimizes the model on the entire data, while CECH and GMMH search for the optimal solutions in each half space individually. Most importantly, since we are not able to directly control the inconsistency level of these methods, they might lead to high inconsistency with initial classification, even if each model is optimized on a separate class (see Forest and Stock data sets). Similar argument holds for classical GMM – although it should allow for optimal fit to the data, it does not take into account the decision boundary between classes.

6.2 Subgroups detection

Decision boundary usually delivers some meaningful information about true structure of clusters. For example, the User Knowledge Modeling data set [13] distinguishes users with very low, low, middle and high knowledge about a given subject. If we knew a coarse separation of the users into two basic classes, {very low, low} and {middle, high}, it should be easier to detect their exact level of knowledge. In this experiment, we want to verify how the information of binary classification influences the clustering results.

To simulate the above scenario, where a decision boundary is closely related with the expected clustering structure, we considered four UCI data sets. For each one we applied the following procedure:

• •

Given $k$ reference groups $Y_{1},\ldots,Y_{k}$ of a data set $X\subset\mathbb{R}^{N}$ we created two classes $Y_{-},Y_{+}$ by merging selected groups together.

• •

We trained SVM classifier on 15% elements drawn randomly from $Y_{-}$ and $Y_{+}$, which induced a linear decision boundary $H$ dividing $X$ into two classes $X_{-}$ and $X_{+}$.

The goal is to discover the reference grouping $Y_{1},\ldots,Y_{k}$. Table 2 contains detailed information about data sets and classes $Y_{-},Y_{+}$.

Given such prepared data sets, we ran all the methods applied in previous experiment. Additionally, we used a version of GMM enhanced with pairwise cannot-link constraints, which is referred as c-GMM [29]. Cannot-link constraints specify the pairs of elements that should not be included into the same group, which suits perfectly to this clustering task555The introduction of must-link constraints is not suitable in this case.. To generate a set of pairwise constraints containing similar knowledge to the decision boundary $H$, we went over all pairs of labeled data points and generated a cannot-link constraint, if one element belonged to $Y_{-}$ and the second belonged to $Y_{+}$.

To compare C3L with other methods we used only two leakage levels $\alpha=0.01$ and $\alpha=0.05$. This choice was motivated by a typical approach used in hypothesis testing, where the significance level is commonly set to $0.01$ or $0.05$. C3L, CEC and CECH were initialized with twice the correct number of clusters (and they were allowed to reduce redundant groups). GMM-based methods were run with the correct numbers of clusters and can thus be expected to perform better than C3L, especially that GMMs describe the clusters by arbitrary Gaussian distributions. The similarity between the obtained clusterings and the ground truth partition $Y_{1},\ldots,Y_{k}$ was evaluated using Normalized Mutual Information (NMI) [3]. NMI is bounded from the above by the value 1, which is attained for identical partitions.

The results presented in Table 3 show that C3L performed better than other methods except the User Knowledge Modeling data set, where GMMH gave very good result. This confirms that working with the entire data set is usually more profitable than finding subgroups in each half space individually (as GMMH and CECH do). Moreover, lower leakage level $\alpha=0.01$ usually led to higher NMI than $\alpha=0.05$. It might follow from the fact that a decision boundary was constructed based on correctly labeled data and, therefore, it was very accurate.

To investigate the influence of the accuracy of decision boundary on the clustering results, we used 15% of data drawn from $Y_{-},Y_{+}$ and assigned incorrect labels to a fixed percentage of them (we considered 0%, 10%, 20% and 30% of erroneous labels). The more labels were misspecified the worse a decision boundary should be.

Table 4 presents the correlation between the accuracy of decision boundary and the normalized mutual information of clustering. One can observe that CECH and GMMH are more sensitive to incorrect decision boundary than both parameterizations of C3L. In consequence, these methods should not be used if there is a risk of unreliable information of class labels. Higher robustness of C3L could be explained by the fact that this model has an access to the entire data set, not only to its part (as GMMH and CECH).

6.3 Improving clusters boundaries

Projection pursuit (PP) is a technique used for analyzing multivariate data by finding its interesting low dimensional views. One dimensional projection is a common choice, which was widely analyzed in the literature [11, 18]. Projections are usually determined by maximizing non-gaussianity. Pena and Prieto [26] showed that minimizing the kurtosis coefficient implies maximizing the bimodality of the projections, which in consequence is useful for detecting clusters. Since clustering in one dimensional subspaces can only generate linear decision boundaries between clusters, PP cannot discover complex data patterns. We will show that given linear separation of data obtained by PP, the use of C3L allows to detect more accurate shapes of clusters.

We considered Statlog data set retrieved from UCI repository concerning credit card applications [15]. Each example is represented by 14 attributes and belongs to one of two classes. First class contains 307 instances, while the second one has 383 objects.

We used R package REPPlab666https://github.com/cran/REPPlab/, which implements several indices for projecting the data on the associated one-dimensional directions. Since we aim at finding clusters, we chose such a direction, which minimizes the kurtosis coefficient of projection. Figure 6(a) presents density estimated from two underlying classes in an optimal one-dimensional subspace. Let us observe that these classes cannot be separated in the reduced space. Given one-dimensional view we applied classical GMM to detect two clusters (see Figure 6). The agreement between true classification and GMM clustering was measured by NMI, which gave the score of 0.24.

GMM clustering in one dimensional space generated linear decision boundary $H$. In order to improve this clustering, we passed $H$ to C3L, which was run with two clusters for four different leakage levels $\alpha\in\{0.01,0.05,0.25,0.5\}$.

The results presented in Figure 6 show that applying C3L with low leakage levels led to the improvement of linear decision boundary produced by GMM in projected space. High values of NMI for $\alpha\in\{0.01,0.05\}$ and its low values for $\alpha\in\{0.25,0.5\}$ prove that the information retrieved by applying PP was meaningful. In other words, strictly unsupervised model-based clustering could not find true structure of classes, while combining one-dimensional projection with constrained clustering gave significantly better results.

6.4 Detection of chemical classes

Finally, we considered a real data set of 2497 chemical compounds which was manually clustered into four categories by the experts in the field [39]. Each compound was characterized by its structural features using Klekota-Roth fingerprint (4860 attributes) [14, 31]. To reduce a dimensionality of the space, PCA was applied to attribute vectors and only five principle components were used.

Additionally, every compound was assigned to active or inactive class based on its binding constant $K_{i}\geq 0$ measured for 5-HT1A receptor, one of the proteins responsible for the regulation of central nervous system: compounds with $K_{i}<50$ were considered active while those with $K_{i}\geq 50$ were treated as inactives [23]. Summarizing, this data set is contained in $\mathbb{R}\times\mathbb{R}^{5}$ space and a decision boundary is given by $H=\{50\}\times\mathbb{R}^{5}$. We investigate whether the information of compounds activity allows to obtain a partition which is more similar to the expert reference grouping with four chemical categories. Observe that this case study represents more realistic scenario than previously prepared experiments.

Figures 7 and 7 show that C3L gave the second best model in terms of BIC (worse than GMM). Moreover, the partition obtained by C3L was the most similar to the reference grouping. Observe that high NMI values coincide with the stabilization of BIC (compare Figure 7 with Figure 7). The highest similarity was achieved for $\alpha=0.15$. C3L with lower leakage levels as well as CECH and GMMH gave worse results, because the binding constant does not reflect exactly the chemical classes (as it was prepared in previous experiment). Since the activity barrier represents a kind of noisy decision boundary it should be taken into account with lower confidence level (higher leakage). This case study demonstrated that C3L can be successfully applied in natural machine learning problems.

6.5 Summary of the experiments

Experimental study showed that C3L can be successfully applied in a wide range of problems. If we have a coarse categorization of data into two basic classes, C3L can be used to construct clusters, which agree with both expert categorization and true data distribution. In other words, C3L is capable of detecting interesting subgroups, which belong to one of two classes with a fixed probability of error. The quality of such model was verified by applying internal clustering measures such as BIC (Section 6.1) as well as by comparing the results with reference grouping created by a domain expert on various examples of data (Section 6.2).

In the case of partial labeling, where only a small sample of data is categorized, we can first apply any binary classifier to construct approximated decision boundary on the entire data space (Section 6.2). If constructed classification is meaningful (accurate), then subgroups can be found by applying C3L with low leakage levels, while in more uncertain situations the leakage levels should be higher. In particular, we showed that compounds activity delivers small amount of information about structural division of chemical space (Section 6.4).

Finally, we demonstrated that C3L can also be used in strictly unsupervised cases. C3L can improve the results of simpler clustering techniques by introducing nonlinearities in clusters descriptions. In particular, we showed that combining C3L with projection pursuit allows to construct better clustering structure than using projection pursuit or model-based clustering individually (Section 6.3).

7 Conclusion and future work

The paper presented a clustering model, C3L, which integrates information coming from the initial classification with the true structure of data. The idea is based on retrieving Gaussian-like clusters which are contained in one of initial classes within a predefined confidence level. Experimental results prove that our algorithm provides high quality model, which can be used for discovering natural subgroups in partially classified data spaces. In the optimization procedure we restricted the problem to special type of Gaussian densities, which is the main limitation of our model. It is worth to eliminate this assumption in future and extend C3L to arbitrary Gaussian distributions using either analytical calculations or by applying some numerical procedures.

We also plan to examine its practical usefulness in real-life problems. In particular, we will focus on applying C3L in text translation systems to find clusters, which combine domain knowledge with a distribution of data (see introduction for details). Since Gaussian mixture model does not fit well to high dimensional data, such as text, it might be needed to extend the model to other probability models or to use projections to lower dimensional spaces. In our opinion, combining C3L with projection pursuit is an approach of great practical potential, which needs further studies. While projection pursuit allows to select optimal low dimensional view of data for finding general clusters regions, C3L is able to use this knowledge to discover detailed clusters description.

Appendix A Algorithm

Given a density model $\mathcal{G}_{1,N-1}$ and a linear constraint $(H,\alpha)$, where $H=\{0\}\times\mathbb{R}^{N-1}$, the one cluster C3L cost function can be evaluated, as a sum of partial costs (see Corollary 1):

[TABLE]

The optimization of $(N-1)$-dimensional density $g_{N-1}\in\mathcal{G}_{N-1}$ is independent of the constraint. Its optimal parameters are the maximum likelihood estimators (MLE) of a mean and covariance of a cluster, i.e.,

[TABLE]

The constraint only affects the form of remaining one dimensional density $g_{1}\in\mathcal{G}_{1}$. Making use of the results and the notations of Theorem 1 it is calculated as:

[TABLE]

The one cluster cost functions of every group are plugged into the expression (8) and determine the overall C3L cost. Its minimization can be performed in an iterative procedure which is a modified online Hartigan algorithm employed in k-means method [10]. Basically, the procedure consists of two steps: initialization and iteration. In the initialization stage, $k\geq 2$ nonempty groups are formed randomly. Then the elements are reassigned between clusters in order to minimize the criterion function. Presented clustering procedure is non-deterministic and one of the local minima is found [12]. To provide more stable and accurate results, the algorithm has to be run a couple of times and a partition with a minimal cost should be chosen. A pseudocode of C3L is given below:

1:  INPUT:

2:  $X\subset\mathbb{R}^{N}$

3:  $k$ - number of clusters

4:  $(H,\alpha)$ - linear constraint

5:  OUTPUT:

6:  Final partition $\mathcal{Y}$ of $X$

7:  INITIALIZATION:

8:  $X=T_{H}(X)$ // data transformation such that $H$//  is mapped into $(\{0\}\times\mathbb{R}^{N-1})$

9:  $\mathcal{Y}\leftarrow$ random partition of $X$ into $k$ groups

10:  for all $Y\in\mathcal{Y}$ do

11:     $g^{Y}_{N-1}\leftarrow MLE(Y{[2:N]}\|\mathcal{G}_{N-1})$ // use standard MLE to find optimal density (19)

12:     $g^{Y}_{1}\leftarrow ConstrMLE(Y{[1]}\|\mathcal{G}_{1}\text{ s.t. }(\{0\}\times\mathbb{R}^{N-1},\alpha))$ // use Theorem 1 for density estimation (20)

13:     $g^{Y}\leftarrow g^{Y}_{1}\cdot g^{Y}_{N-1}$

14:  end for

15:  ITERATION:

16:  while NOT Done do

17:     $Done\leftarrow True$

18:     for all $x\in X$ do

19:        $Y_{new}\leftarrow arg\max\limits_{Y\in\mathcal{Y}}\Delta E_{\{0\}\times\mathbb{R}^{N-1}}^{\alpha}(x,Y)$ // find a membership of $x$ maximizing the decrease of cost

20:        if $Y_{new}\neq x.cluster$ then

21:           $Done\leftarrow False$

22:           $Reassign(x,Y_{new},Y_{old})$ // reassign $x$ from $Y_{old}$ to $Y_{new}$

23:           $Update((Y_{new},g^{new}),(Y_{old},g^{old}),x,\alpha)$ // recalculate clusters parameters after the change

24:        end if

25:     end for

26:  end while

Appendix B Proof of Lemma 1

We consider the limiting case of $\alpha\to 0$. Therefore, without loss of generality we may assume that there exists $\alpha_{0}>0$ such $m_{X}^{\alpha}$ and $\sigma_{X}^{\alpha}$, for $\alpha<\alpha_{0}$, are given by (15), i.e.,

[TABLE]

where

[TABLE]

Let us calculate the limiting value of $m_{X}^{\alpha}$:

[TABLE]

Making use of Taylor expansion we have

[TABLE]

where $\xi_{x}\in(0,x)$. Consequently, in our situation there exists $\xi_{\alpha}\in(0,\frac{4}{(p^{\alpha})^{2}}(1+\frac{\hat{\sigma}_{X}^{2}}{\hat{m}_{X}^{2}}))$ such that for $\varepsilon_{\alpha}=-\frac{1}{8}(1+\xi_{\alpha})^{-\frac{3}{2}}\xi_{\alpha}^{2}$ we get

[TABLE]

Clearly, (I) converges to $2(\hat{m}_{X}+\frac{\hat{\sigma}_{X}^{2}}{\hat{m}_{X}})$, as $\alpha\to 0$. We consider the term (II). Observe that for sufficiently small $\alpha>0$, we have:

[TABLE]

Concluding $m_{X}^{\alpha}\to\hat{m}_{X}+\frac{\hat{\sigma}_{X}^{2}}{\hat{m}_{X}}$, as $\alpha\to 0$.

From the above calculation we directly get,

[TABLE]

which completes the proof.

Appendix C Proof of Lemma 2

Let $\alpha>0$. The violation of the linear constraint $(\{0\},\alpha)$ by a density $\mathcal{N}(m,\sigma)$ is verified by the condition:

[TABLE]

Theorem 1 states that the parameters of optimal one dimensional Gaussian density are calculated as MLE within the cluster if only it does not violate the linear constraint. However, in the limiting case there exists $\alpha_{0}$ such that for all $\alpha<\alpha_{0}$ the constraint is violated and consequently the formulas (15) for $m_{\pm}^{\alpha},\sigma_{\pm}^{\alpha}$ give optimal solutions.

Let $x>0$. For normal distributions $g_{\pm}^{\alpha}=\mathcal{N}(m_{\pm}^{\alpha},\sigma_{\pm}^{\alpha})$ we have

[TABLE]

Making use of equality $\sigma_{\pm}^{2}=\frac{m_{\pm}^{2}}{(p^{\alpha})^{2}}$ we get

[TABLE]

and consequently

[TABLE]

Since $m_{\pm}^{\alpha}\to\hat{m}_{\pm}+\frac{\hat{\sigma}_{\pm}^{2}}{\hat{m}_{\pm}}$, as $\alpha\to\infty$ (see Theorem 2) then

[TABLE]

Therefore,

[TABLE]

iff

[TABLE]

The inequality (21) holds iff

[TABLE]

In the limiting case the last inequality expands to:

[TABLE]

which completes the proof.

Acknowledgement

This research was partially supported by the National Science Centre (Poland) grant no. 2016/21/D/ST6/00980 and grant no. 2015/19/B/ST6/01819.