Keywords: classication non supervisée, contraintes utilisateur, Programmation par Contraintes, clustering bi-critère constrained clustering, bicriterion clustering, constraint programming, modelling constrained clustering, bicriterion clustering, constraint programming, modelling, ltering algorithm

La classication non supervisée, souvent appelée par le terme anglais de clustering, est une tâche importante en Fouille de Données, pour lesquelles de nombreuses approches ont été développées. Depuis une dizaine d'années, la classication non supervisée a été étendue pour intégrer des contraintes utilisateur permettant de modéliser des connaissances préalables dans le processus de clustering. Diérents types de contraintes utilisateur peuvent être considérés, des contraintes pouvant porter soit sur les clusters, soit sur les instances. Par exemple, des contraintes peuvent exprimer des bornes sur la taille ou le diamètre des clusters. Des contraintes sur des instances expriment que deux objets doivent ou ne doivent pas être dans le même cluster (contraintes must-link et cannot-link).

Dans cette thèse, nous étudions le cadre de la Programmation par Contraintes (PPC) pour modéliser les tâches de clustering sous contraintes utilisateur. Des avancés récentes en PPC ont rendu ce paradigme puissant pour résoudre des problèmes combinatoires. Les principes de la PPC sont : (1) le programmeur spécie le problème d'une façon déclarative par un Problème de Satisfaction de Contraintes ; (2) le solveur recherche la (les) solution(s) en intégrant la propagation de contraintes à l'exploration de l'espace de recherche. Utiliser la PPC a deux avantages principaux : la déclarativité, qui permet d'intégrer aisément des contraintes utilisateur et la capacité de trouver une solution optimale qui satisfait toutes les contraintes (s'il en existe).

Nous proposons deux modèles basés sur la PPC pour le clustering sous contraintes utilisateur. Les modèles sont généraux et exibles, ils permettent d'intégrer des contraintes d'instances must-link et cannot-link et diérents types de contraintes sur les clusters. Ils orent également à l'utilisateur le choix entre diérents critères d'optimisation. An d'améliorer l'ecacité, divers aspects sont étudiés. Les expérimentations sur des bases de données classiques et variées montrent qu'ils sont compétitifs par rapport aux approches exactes existantes. Nous montrons que nos modèles peuvent être intégrés dans une procédure plus générale et nous l'illustrons par la recherche de la frontière de Pareto dans un problème de clustering bi-critère sous contraintes utilisateur. A notre meilleure connaissance, notre approche est la première pour les problèmes de classication non supervisée sous contraintes utilisateur, qui intègrent le critère de diamètre des clusters, de la marge entre clusters, de la somme des dissimilarités intra-cluster et du bi-critère marge-diamètre, en présence de contrainte utilisateur.

Introduction

Ce chapitre synthétise le contexte et la problématique de la thèse. Nous présentons brièvement les notions principales dans le mémoire: la classication non supervisée (clustering), le clustering sous contraintes utilisateur et la Programmation par Contraintes (PPC). Dans cette thèse, nous étudions le cadre de la PPC pour modéliser les tâches de clustering sous contraintes utilisateur. Nous proposons deux modèles basés sur la PPC pour le clustering sous contraintes utilisateur. Les modèles sont généraux et exibles, ils permettent d'intégrer des contraintes d'instances must-link et cannot-link et diérents types de contraintes sur les clusters. Ils orent également à l'utilisateur le choix entre diérents critères d'optimisation. An d'améliorer l'ecacité, divers aspects sont étudiés. Le chapitre précise les contributions, les publications et l'organisation du mémoire. Cluster analysis, also called clustering, is a fundamental task in Data Mining, which aims at partitioning a given set of objects into homogeneous and well-separated subsets called classes or clusters. Homogeneity means that the objects inside the same cluster must be similar while well separation expresses that objects in a cluster must be dierent from the objects belonging to other clusters. The problem has been studied since the 18 th century with applications in natural sciences, economics, sociology, geology and many other elds. Cluster analysis is an example of unsupervised learning, as it does not require predened category labels of objects. It helps to nd and understand the nature or the structure of the data. For instance, in marketing, by analysing the buying records of supermarkets, cluster analysis can nd groups of customers with similar behaviour. Cluster analysis is a dicult task with many challenges, for instance: how to dene the similarities between objects? how to decide the number of clusters? or how to evaluate the clusters? Although there have been hundreds of dierent approaches in the literature, such a general clustering algorithm does not exist. In cluster analysis, there is usually no single correct result and a clustering algorithm may give good results with some types of data and bad results with other data. In many real applications, there exists some background knowledge about the data that could be useful in clustering. Since the last decade, cluster analysis has been extended to semi-supervised clustering, so as to integrate previous knowledge on data to clustering algorithms. The background knowledge is usually given in the form of user constraints. The constraints can be classied into cluster-level constraints, specifying requirements on the clusters, or instance-level constraints, specifying requirements on pairs of objects. Most of the attention has been put on instance-level constraints, also called pairwise constraints, rst introduced in [START_REF] Wagsta | [END_REF]. The instance-level constraints are composed of must-link and cannot-link constraints. A must-link constraint on two data objects species that they have to appear in the same cluster, whereas a cannot-link constraint species that the two objects must not be in the same cluster. Over the last decade, many works have been done to extend classical algorithms for handling must-link and cannot-link constraints, as for instance extensions of COBWEB [START_REF] Wagsta | [END_REF], k-means [START_REF] Wagsta | [END_REF], Bilenko et al. 2004], hierarchical clustering [Davidson & Ravi 2005a] or spectral clustering [START_REF] Lu & Carreira-Perpinan | [END_REF][START_REF] Wang | [END_REF], etc. This is achieved either by modifying the dissimilarity measure, or the objective function, or the search strategy. However, there is no general solution to extend traditional algorithms to dierent types of constraints.

Recently, there has been a growing interest in developing declarative framework for Data Mining [START_REF] Raedt | [END_REF], including the problem of constrained clustering.

Declarative frameworks may be less ecient than classic algorithms for a specic task, but it is more exible and general, easy to incorporate new knowledge to the task. [Guns et al. 2013] proposes a framework in Constraint Programming for k-pattern set mining. [START_REF] Métivier | [END_REF]] presents a constraint-based language for expressing queries to discover patterns in Data Mining. A Satisability of boolean formula (SAT) framework for constrained clustering is proposed in [Davidson et al. 2010] for cluster analysis with 2 clusters. The framework integrates dierent kinds of user-constraints: must-link, cannotlink, diameter and split constraints. [Mueller & Kramer 2010] proposes an approach to constrained clustering based on Integer Linear Programming. This approach takes a set of candidate clusters as input and constructs a clustering by selecting a suitable subset. It allows dierent kinds of constraints on clusters or on the set of clusters, but no constraint on individual objects. [START_REF] Babaki | [END_REF] proposes an exact approach for constrained clustering with the criterion of minimizing the within-cluster sum of squares, based on Integer Linear Programming. It allows must-link, cannot-link and all constraints that are anti-monotone.

Among generic optimization tools, Constraint Programming (CP) is a powerful paradigm for solving combinatorial search problems that relies on a wide range of techniques from Articial Intelligence, Computer Science and Operational Research. The basic idea of Constraint Programming is to solve problems by stating constraints that must be satised by the solution. The main principles in CP are: (1) users specify declaratively the problem as a Constraint Satisfaction Problem; (2) solvers search for solutions by constraint propagation and search. We claim that Constraint Programming might be a useful tool for the task of cluster analysis with constraints. Relying on Constraint Programming has two main advantages: the declarativity, which enables to easily add new constraints and the ability to nd an optimal solution satisfying all the constraints (when there exists one). In this dissertation, we present our research on Constraint Programming to solve the problem of Constrained Clustering.

1.1. Dissertation Contributions 3

Dissertation Contributions

The main contributions of this dissertation are three folds:

• We propose two models based on Constraint Programming, allowing to nd an optimal solution for clustering under constraints, given an optimization criterion.

• We show that such a framework can easily be embedded in a more general process and we illustrate this on the problem of nding the Pareto front of a bi-criterion (split-diameter) optimization process.

• We have proposed new global constraints to improve the performance of the framework.

We presented the rst model [Dao et al. 2013d, Dao et al. 2013c, Dao et al. 2013a],

based on Constraint Programming, which enables to design clustering tasks by specifying an optimization criterion and some constraints either on the clusters or on pairs of objects. We proposed improvements of the search strategy and a ltering algorithm [Dao et al. 2013b] in order to enhance the performance of the model. In our framework, several classical optimization criteria are considered and they can be coupled with different kinds of constraints. Experimental results show that our model dominates the state-of-the-art exact clustering algorithms for the criterion of diameter.

Next, we analysed the limitation of the framework and we proposed a new and improved model in Constraint Programming [Dao et al. 2014a, Dao et al. 2014c, Dao et al. 2014b]. The second model is more exible, as it does not require to set the number of clusters. Moreover, the model is lighter in terms of the number of variables and constraints. We also developed dedicated global constraints for two criteria. Experiment results show that the second model it not only more general but also much more ecient. Finally, we studied the problem of constrained clustering with a multi-criterion optimization and demonstrated the use of our model in order to solve the problem of clustering with the bi-criterion of diameter and split [Dao et al. 2014a, Dao et al. 2014c, Dao et al. 2014b].

List of publications

Dissertation Organization

This dissertation is organized as follows:

Chapter 2 is dedicated to preliminaries on Constraint Programming.

Chapter 3 gives background on the problem of cluster analysis, mainly from an optimization viewpoint. We present dierent approaches of cluster analysis that are related to our research. Finally, we present the state of the art of constrained clustering.

Chapter 4 brings an overview of recent declarative frameworks for Data Mining, specially those that are related to cluster analysis.

Chapter 5 presents our rst model for solving the problem of constrained clustering.

The model, based on CP, supports instance-level constraints, specifying requirements on pairs of objects, as well as dierent cluster-level constraints, specifying requirements on clusters. It also allows users to choose among dierent optimization criteria. The contributions of this Chapter have been published in [Dao et al. 2013d Chapter 7 gives the introduction of multi-objective optimization, and demonstrate the use of our model for nding the optimal Pareto solutions of a bi-criterion optimization process. The contributions of this Chapter are in [Dao et al. 2014a, Dao et al. 2014c, Dao et al. 2014b].

In the conclusion the dissertation is summarized and perspectives are discussed.

Chapter 2

Constraint Programming

Dans ce chapitre, nous présentons des connaissances de base en Programmation par Contraintes (PPC): le Problème de Satisfaction de Contraintes, la propagation de contraintes, les contraintes globales et le ltrage, la recherche de solutions et la notion de symétrie. Un Problème de Satisfaction de Contraintes (CSP) est un triplet X, Dom, C où X = x 1 , x 2 , . . . , x n est un n-tuple de variables, Dom = Dom(x 1 ), Dom(x 2 ), . . . , Dom(x n ) est un n-tuple de domaines (x i ∈ Dom(x i )), C = C 1 , C 2 , ..., C t est un t-tuple de contraintes où chaque contrainte C i exprime une condition sur un sous-ensemble de X. Une solution d'un CSP est une aectation complète de valeurs a i ∈ Dom(x i ) à chaque variable x i qui satisfait toutes les contraintes de C. La propagation de contraintes opère sur une contrainte c et enlève du domaine des variables de c des valeurs dont on est certain qu'elles ne pourront pas apparaître dans une solution. Les contraintes globales représentent des relations sur des ensembles de variables et en général, elles peuvent être exprimées par une conjonction de contraintes élémentaires. Cependant, dénies d'une façon globale, les contraintes globales bénécient de mécanismes de propagations beaucoup plus ecaces, qui sont eectuées par des algorithmes de ltrage. Un solveur de contraintes itère une étape de propagation suivie d'une étape de branchement jusqu'à ce qu'une solution soit trouvée. This chapter gives an introduction and basic notions of Constraint Programming (CP).

A large part of this chapter is based on the books [Apt 2003] and [START_REF] Rossi | [END_REF]].

General Concepts

In this section, we present general concepts related to Constraint Programming. Constraint Programming is a powerful paradigm for solving combinatorial search problems that relies on a wide range of techniques from Articial Intelligence, Computer Science and Operational Research. The main principles in CP are: (1) users specify declaratively the problem in a Constraint Satisfaction Problem; (2) solvers search for solutions by constraint propagation and search.

Denition 1 Let x be a variable. The domain of x, denoted by Dom(x), is a set of possible values that can be assigned to x. For each variable x, only a single value in Dom(x) can be assigned to x.

Example 1 x 1 , x 2 , x 3 are variables with the domains: Dom(x 1 ) = {1, 2, 3}, Dom(x 2 ) = {2, 5} and Dom(x 3 ) = {1, 3, 5}.

Denition 2 Let X = {x 1 , x 2 , ..., x n } be a set of variables. An assignment is a n-tuple t = (t 1 , . . . , t n ) ∈ Dom(x 1 ) × . . . × Dom(x n ) such that x i = t i for all x i ∈ X, n = |X|.

Example 2 The tuple (1, 2, 5) is an assignment of variables given in Example 1. It is equivalent with x 1 = 1, x 2 = 2 and x 3 = 5.

Denition 3 Let X = {x 1 , x 2 , ..., x n } be a set of variables. A constraint C on an ordered tuple of variables, denoted by V ar(C) = (x i1 , . . . , x ik ) ⊆ X, is a subset of the cartesian product of the domains of the variables in V ar(C), i.e., C ⊆ Dom(x i1 ) × . . . × Dom(x ik ), that species the allowed combinations of values for the variables. V ar(C) is the scope of the constraint and |V ar(C)| is the arity of the constraint C. C is called an unary constraint if it is dened on one variable and is called a binary constraint if it is dened on two variables. In many cases, constraints can be expressed by mathematical formulas on variables.

Example 3 Given the variables dened in Example 1, we consider two binary constraints: C 1 ≡ x 1 = x 2 and C 2 ≡ x 1 + x 3 ≤ 6. The rst constraint can be dened as {(1, 2), (1, 5), (2, 5), (3, 2), (3, 5)}. The second constraint is equivalent to {(1, 1), (1, 3), (1, 5), (2, 1), (2, 3), (3, 1), (3, 3)}.

Denition 4 A Constraint Satisfaction Problem (CSP) P is a triple X, Dom, C where:

• X = x 1 , x 2 , . . . , x n is a n-tuple of variables,

• Dom = Dom(x 1 ), Dom(x 2 ), . . . , Dom(x n ) is a corresponding n-tuple of domains such that x i ∈ Dom(x i ),

• C = C 1 , C 2 , ..., C t is a t-tuple of constraints where each constraint C i expresses a condition on a subset of X.

Example 4 Based on Examples 1 and 3, we have a CSP dened as:

x 1 ∈ {1, 2, 3}, x 2 ∈ {2, 5}, x 3 ∈ {1, 3, 5}

x 1 = x 2 , x 1 + x 3 ≤ 6

2.1. General Concepts 9 Denition 5 Given a CSP P = X, Dom, C with X = {x 1 , x 2 , ..., x n }, Dom = Dom(x 1 ), Dom(x 2 ), . . . , Dom(x n ) and a set of constraints C = {C 1 , C 2 , ..., C m }, a tuple t = (t 1 , . . . , t n ) satises a constraint C i dened on the variables x i 1 , . . . , x i k if the tuple (t i 1 , . . . , t i k ) ∈ C i .

Example 5 Given the CSP in Example 4, the tuple (1, 2, 3) satises the constraint x 1 = x 2 since the tuple (1,2) (which is equivalent to the assignment x 1 = 1, x 2 = 2) satises the constraint x 1 = x 2 . Denition 6 Given a CSP P = X, Dom, C with X = {x 1 , x 2 , ..., x n }, Dom = Dom(x 1 ), Dom(x 2 ), . . . , Dom(x n ) and a set of constraints C = {C 1 , C 2 , ..., C m }, a tuple t = (t 1 , . . . , t n ) ∈ Dom(x 1 ) × . . . × Dom(x n ) is a solution of this CSP if it satises every constraint in C. Example 6 Given the CSP in Example 4, a solution to the problem is the tuple (1,2,3).

The assignment of x 1 = 1, x 2 = 2 and x 3 = 3 satises all the constraints in this CSP.

Denition 7 A constraint optimization problem (COP) is a CSP P = X, Dom, C together with an objective function f : Dom(x 1 ) × . . . × Dom(x n ) → R to be optimized. An optimal solution to a COP is a solution to P that optimizes the function f . Example 7 Given the CSP in Example 4 and a function f = x 1 +x 2 +x 3 to be maximized, we have a COP dened as:

x 1 ∈ {1, 2, 3}, x 2 ∈ {2, 5}, x 3 ∈ {1, 3, 5}

x 1 = x 2 , x 1 + x 3 ≤ 6 maximize x 1 + x 2 + x 3
An optimal solution of this COP is the tuple (1,5,5), which is equivalent to the assignment x 1 = 1, x 2 = 5, x 3 = 5. Another optimal solution of this COP is the tuple (3,5,3).

In general, CSPs are NP-hard. However, techniques used in CP solvers allow to solve eciently many dicult problems. For solving a practical problem by CP, the problem must be modeled as a CSP or a COP, which means it is needed to precise:

• Denition of variables and theirs domains.

• Specication of the constraints between variables.

• Denition of the objective function on related variables if it is a COP.

CSPs and COPs can be modeled and then solved by a constraint solver such as: Gecode [Gecode Team ], Choco [choco Team 2010], Comet [ Van Hentenryck & Michel 2005]

Constraint Propagation

In general, a constraint solver nds one or all solutions of a CSP by using successively constraint propagation and search. Constraint propagation consists of removing inconsistent values from domains of variables. In consequence, the search space can be signicantly reduced.

Denition 8 Let C be a constraint dened on variables V ar(C) = {x i1 , . . . , x im }. A support for a value v ∈ Dom(x ij ) is an assignment t = (t i1 , . . . , t im ) such that t ij = v and (t i1 , . . . , t im ) ∈ C. If there doesn't exist such a support, the value v is inconsistent with the constraint C.

The constraint propagation phase is proceeded by propagation algorithms, or also called ltering algorithms.

Denition 9 A ltering algorithm associated with a constraint C is an algorithm that removes inconsistent values of variables involved in the constraint. If the ltering algorithm cannot remove any inconsistent value in the domain of variables, the constraint C is locally consistent.

A constraint is locally consistent if the ltering algorithm cannot remove any inconsistent values. In constraint propagation, ltering algorithms are applied repeatedly until all constraints are locally consistent, with respect to domains of variables. When this process terminates, the CSP is called locally consistent. Many types of local consistencies have been studied in Constraint Programming, including node consistency, arc consistency, generalized arc consistency and path consistency. The simplest type is the node consistency which requires that every unary constraint C dened on a variable x i is satised by all values in Dom(x i ). In practice, the node consistency can be achieved by reducing the domain of each variable to the values that satisfy all unary constraints. As a result, unary constraints can be handled one time at the beginning and be neglected later.

Denition 10 A CSP is node consistent if for every variable x ∈ X, every unary constraint dened on x is consistent with the domain of x.

Example 8 Given the CSP in Example 4 with a new unary constraint

C 3 ≡ x 1 ≤ 2.
The value 3 is inconsistent in Dom(x 1 ) with the constraint C 3 . By ltering this value, Dom(x 1 ) = {1, 2} and the CSP is node consistent.

The next level of consistency is arc consistency, which is dened on binary constraints.

Denition 11 Given a binary constraint C dened on variables x and y, with the domains Dom(x) and Dom(y) respectively. The constraint C is arc consistent (AC) if for every value a ∈ Dom(x), there exists a value b ∈ Dom(y) such that the assignment x = a, y = b satises the constraint C and vice versa:

∀a ∈ Dom(x), ∃b ∈ Dom(y) : (a, b) ∈ C ∀b ∈ Dom(y), ∃a ∈ Dom(x) : (a, b) ∈ C
A CSP is arc consistent if all its binary constraints are arc consistent.

Constraint Propagation 11

Example 9 The two binary constraints of the CSP given in Example 4 are arc consistent.

Considering a new binary constraint: C 3 ≡ x 1 = x 3 , it is not arc consistent because if we assign the value 1 to the variable x 1 , there is no value for x 3 in its domain to satisfy the constraint x 1 = x 3 .

The notion of arc consistency can be generalized to arbitrary constraints.

Denition 12 A constraint C is generalized arc consistent (GAC) if for each variable involved in the constraint, for any of its value there exists a value for every other variables such that the assignment of those values to variables satises the constraint.

In another way, a constraint C is generalized arc consistent if every value of every variable involved in the constraint has at least a support. A CSP is generalized arc consistent if all of its constraints are generalized arc consistent.

The basic arc consistency algorithm is presented in Algorithm 1. All the binary constraints are revised each time an inconsistent value is detected and removed. In the literature, many other arc consistency algorithms (AC-2, . . . , AC-7) have been proposed to reduce the complexity.

Algorithm 1: Basic Arc Consistency algorithm 1 Function evise@x i , x j A 2 DELETED ← false 3 foreach value a i in Dom(x i ) do 4 if there is no a j in Dom(x j ) such that (a i , a j ) satises all binary constraints on (x i , x j ) then Example 10 Assume the variables x 1 , x 2 and x 3 taking integer values in [1, 4]1 with constraints: x 1 < x 2 and x 2 < x 3 . Suppose that arcs will be revised in the order (x 1 , x 2 ), (x 2 , x 1 ), (x 2 , x 3 ), (x 3 , x 2 ).

• we rst select the arc (x 1 , x 2 ) and remove the value 4 from Dom(x 1 ) because if

x 1 = 4, there exists no value for x 2 such that x 1 < x 2

• arc (x 2 , x 1 ) is revised and the value 1 is removed from Dom(x 2 ) for the same reason.

• arc (x 2 , x 3 ) is revised and the value 4 is removed from Dom(x 2 ). Now Dom(

x 1 ) = [1, 3], Dom(x 2 ) = [2, 3], Dom(x 3 ) = [1, 4]
• arc (x 3 , x 2 ) is revised and the values 1 and 2 are removed from Dom(x 3 ). Now

Dom(x 1 ) = [1, 3], Dom(D 2 ) = [2, 3], Dom(x 3 ) = [3, 4]
• all the arcs are revised again: when revising (x 1 , x 2 ), the value 3 is removed from Dom(x 1 )

• Finally Dom(x 1 ) = [1, 2], Dom(x 2 ) = [2, 3], Dom(x 3 ) = [3, 4]
The Arc Consistency techniques reduce signicantly the search space in many CSPs but there still exist many other possible inconsistencies.

Example 11 Assume the variables x 1 , x 2 and x 3 taking integer values in {1, 2} with constraints: x 1 = x 2 , x 2 = x 3 and x 3 = x 1 . The Arc Consistency techniques can not detect any inconsistencies in this case although there is no solution.

Path Consistency is stronger than Arc Consistency, it detects inconsistencies in every path in the constraint graph. It removes more inconsistencies but it is too costly and is not ecient in most CSPs. Moreover, the Path Consistency is not able to detect all of the inconsistencies.

Global constraints

A global constraint is a constraint that expresses a relationship between a non-xed number of variables. This relationship can be expressed by a conjunction of simpler constraints.

Denition 13 Let C = {C 1 , . . . , C m } be a set of constraints. The constraint C G , which is equal to the conjunction of all the constraints of C:

C G = ∧{C 1 , . . . , C m }, is a global constraint.
A typical example of a global constraint is the AllDierent(x 1 , x 2 , . . . , x n ) constraint. This constraint requires that values assigned to the variables x 1 , . . . , x n must be pairwise distinct.

Denition 14 Let X = {x 1 , . . . , x n } be a set of variables. The AllDierent constraint on X is dened as:

AllDierent(x 1 , . . . , x n ) = {(d 1 , . . . , d n ) | ∀ i : d i ∈ Dom(x i ), ∀ i =j : d i = d j }
For instance, the constraint AllDierent(x 1 , x 2 , . . . , x n ) can be replaced by n(n-1) 2 binary constraints:

x 1 = x 2 , x 1 = x 3 , . . . , x n-1 = x n .
Using global constraints makes a CSP become more intuitive. But more importantly, powerful ltering algorithms can be designed because the set of simple constraints can be considered as a whole at the same time. Many ltering algorithms for global constraints Denition 15 Given a graph G = (V, E), a matching in G is a set M ⊆ E of disjoint edges, i.e no two edges in M share a vertex. A matching M is a maximum matching if it contains the largest possible number of edges.

Denition 16 Let X = x 1 , . . . , x n be a sequence of variables with respective domains Dom(x 1 ), . . . , Dom(x n ). The bipartite graph

G = (V, E), with V = X ∪ Dom(x 1 ) ∪ . . . ∪ Dom(x n ) and E = {(x i , d)|x i ∈ X, d ∈ Dom(x i )} is called the value graph of X.
In the value graph G, each variable x i corresponds to a variable-node while each value v corresponds to a value-node. There is an edge between a variable-node x i and a value-node v i the value v is in the domain of x i .

Example 12 Assume we have 3 variables x 1 , x 2 and x 3 with the domains: x 1 ∈ {1, 3}, x 2 ∈ {1, 3}, x 3 ∈ {1, 2, 3}. The value graph of those variables and a maximum matching are shown in Figure 2.1.

Let X = {x 1 , . . . , x n } be a set of variables and let G be the value graph of X, [Régin 1994] introduced two important propositions on the relationship between assignments of X and matchings.

• There is a matching of cardinality n if and only if the AllDierent constraint is satisable.

• An edge (x i , v) belongs to a matching of cardinality n if and only if the value v is consistent with the constraint.

From that observation, the ltering algorithm in [Régin 1994] consists of two steps:

• Determine if there exists at least one matching of cardinality |X|. If such a matching does not exist, the constraint cannot be satised.

• Mark all edges (x i , v) that could not belong to any matching of cardinality n. For each edge (x i , v) satisfying the condition, the value v is inconsistent and can be removed from the domain of x i . The rst step can be solved by nding a matching of maximum cardinality M in the value graph. Because there is only edges between variable-nodes and value-nodes, the maximum cardinality of matchings in the value graph is at most |X|.

For the second step, [Régin 1994] constructs a residual graph G M that is a directed version of the graph G, with the direction based on M . Edges that belong to the maximum matching M are oriented from the value-nodes to variable-nodes. Otherwise, edges are oriented from value-nodes to variable-nodes. The ltering algorithm nds all strongly connected components in G M and marks all edges in those components. A simple depthrst search marks all edges that lie on an even-length path starting at the free nodes. It is proven that those marked edges belong to at least one maximum cardinality matching and the unmarked edges do not belong to the matching M and cannot be a part of any other maximum cardinality matching. Hence, the AllDierent constraint is made generalized arc consistent by removing all the values v ∈ D(x i ) such that (x i , v) is an unmarked edge.

Example 13 Given the value graph in Example 12, the graph G M and marked edges are presented in Figure 2.2. There are two unmarked edges (x 3 , 1) and (x 3 , 3). Those edges do not belong to any maximum cardinality matching, therefore we can remove value 1 and 3 from the domain of x 3 .

Search

Denition 17 Given a CSP P 0 , a search tree for P 0 is a rooted tree such that:

• its nodes are CSPs,

• its root is P 0 ,

• if P 1 , . . . , P m where m > 0 are all direct descendants of P 0 , then the union of the solution sets of P 1 , . . . , P m is equal to the solution set of P 0 .

A node in the search tree is a dead end if it does not lead to a solution. The general idea is that the search tree is built by branching (also called splitting) a CSP into smaller CSPs until each leaf is a dead end or a solved CSP. At each node, the constraint propagation is applied to the CSP until the domains of variables are consistent. There are three cases: • If any of the variable domains is empty, the CSP has no solution and the search explorers other branches.

• Otherwise, the current CSP is locally consistent and it is split into smaller CSPs either by splitting a constraint or by splitting the domain of a variable. In this thesis, we only consider the latter technique.

The strategy of extending a node in the search tree is often called a branching strategy. An example of a search tree is given in Figure 2.3. Solving a CSP can be complete or incomplete:

• A complete search guarantees to nd all existing solutions, or proves that the CSP does not have any solution. A complete search also guarantees to nd the optimal solution of a COP. The chronological backtracking search is the most well known complete approach. A search is completed if any leaf in the search tree is either a dead end or a solved CSP.

• An incomplete search does not explore the search space completely, but focuses on branches which have a high probability of containing a solution. It is often eective at nding a solution of CSPs if one exists. In COPs, an incomplete search can be used to nd an approximation to an optimal solution. Incomplete searches can not be used to prove that the CSP does not have a solution or to prove the optimality of a solution as the searches does not explore all possible branches of the search tree.

Backtracking search in CP

The naive backtracking algorithm begins with the root node P 0 which is the original CSP.

The node is then split following a branching strategy. For example, in the branching strategy of enumeration, the node P 0 is extended by selecting a variable, for example the variable x, and splitting its domain into singleton sets. For each value a i in the domain of x, a descendant node P i is generated. P i is a copy of P 0 except the domain of x is replaced by a singleton: Dom(x) = {a i }. The backtracking algorithm then visits a node P i and the constraint propagation is applied to P i . There are three cases after that:

Chapter 2. Constraint Programming

• The domain of a variable is reduced to an empty set: P i is a failed CSP, the search backtracks to visit other non-visited nodes.

• All the variable domains are singletons: these singleton domains form a solution, the search backtracks to visit other non-visited nodes.

• Otherwise: Node P i is extended and the search visits one of the descendent nodes.

The search is terminated after visiting every generated nodes. Three popular branching strategies often used are:

• Enumeration: A branch is generated for each value in the domain of the variable, as explained above.

• Binary choice points: A variable x and a value a ∈ Dom(x) are chosen. There are two descendant nodes P 1 and P 2 . In P 1 , the domain of x is replaced by a singleton Dom(x) = {a} whereas in P 2 , the domain of x contains all remaining values.

• Domain splitting: Here the domain of the variable is split to create the branching. Suppose that a variable x and a value a is chosen for branching. There are two descendant nodes P 1 and P 2 . In P 1 , the domain of x contains all values which are less or equal to a, and P 2 contains the remaining values.

Strategies for choosing variables and for choosing values for branching are extremely important, since they can reduce drastically the search tree.

Variable and value ordering heuristics

For branching a node, we have to decide how to select a variable and how to split its domain. The order in which variables and values are selected is extremely important in the search process since correct choices can reduce drastically the number of nodes in the search tree. The variable ordering may be either static or dynamic. In the static ordering, the order of the variables is specied before the search begins whereas in the dynamic ordering, the choice of the variable is considered as the search processes. [Haralick & Elliott 1979] suggests that "the best search order is the one which minimizes the expected length or depth of each branch". In order to do that, [Haralick & Elliott 1979] introduced the Fail-First Principle that tries to select the variable which can quickly leads to a dead end node. The underlying idea is to detect failures early to avoid useless branches. If the current node P cannot lead to a solution, then the sooner we detect the failure the better. The most common variable ordering heuristic with the Fail-First Principle is the Smallest Remaining Domain method. In this method, the variable with the smallest domain size is selected for branching with the idea that the variable is most likely to lead to a deadend. This method is based on the assumption that each value in a domain has the same probability to be in a solution, so the variable with the smallest remaining domain size has the smallest chance to be successful. Other heuristics may be used during the search to determine the next variable to be instantiated:

• the variable that participates in a higher number of constraints,

• the variable that has the largest number of constraints with instantiated variables. The order of values for each variable is important too, a correct choice of values can help to reduce the search tree. The example below demonstrates the importance of the variable ordering.

Example 14 Find distinct digits for the letters S, E, N, D, M, O, R, Y such that

S E N D + M O R E = M O N E Y
To model this problem, we can dene a variable with the domain [0, 9] for each letter:

X = {x S , x E , x N , x D , x M , x O , x R , x Y }.
The following constraints are added:

• AllDif f erent(x S , x E , x N , x D , x M , x O , x R , x Y ) • x S = 0, x M = 0 • 1000x S +100x E +10x N +x D +1000x M +100x O +10x R +x E = 10000x M +1000x O + 100x N + 10x E + x Y
By using arc consistency techniques, we can reduce the domain of variables to:

Dom(x S ) = {9}, Dom(x E ) = [4, 7], Dom(x N ) = [5, 8], Dom(x D ) = [2, 8], Dom(x M ) = {1}, Dom(x O ) = {0}, Dom(x R ) = [2, 8], Dom(x Y ) = [2, 8]
Figures 2.4 and 2.5 show two search trees with dierent strategies for the next variable to be instantiated. These search trees are generated by Gist environment of the Gecode solver where a blue circle denotes a stable state but not yet a solution, a red square represents a fail state (there is no solution) and a green diamond is a solution.

Figure 2.5: Search tree with the strategy: the variable with the biggest domain is selected for instantiation

Search for COPs

To solve a COP, a common approach is to use the branch-and-bound technique. A backtracking search is rst used to nd a solution which satises all the constraints without considering the optimization criterion. When a solution is found, the value of the optimization function is computed and an additional constraint is added to the CSP in order to forbid solutions that are not better than this solution. The solver continues to search for solutions of the new CSP and every time a solution is found, an additional constraint is added to the CSP. When the search space is explored completely, the last solution found is proven optimal. Let us consider the following example for illustrating a Constraint Optimization Problem and search strategies.

Example 15 Find an assignment of digits to letters such that

S E N D + M O S T = M O N E Y
and the value of M ON EY is maximized. The optimal solution of this problem is SEN D = 9782, M OST = 1094 and the value M ON EY = 10876. This problem can be modeled as a Constraint Optimization Problem, using eight variables x S , x E , x N , x D , x M , x O , x T , x Y , whose domain is the set of digits [0, 9], and a variable V of domain integer, which represents the objective function, which is to be maximized. Constraints that specify the problem are:

• the letters S and M must not be equal to 0: x S = 0, x M = 0 Figure 2.6: Search trees (blue circle: a stable state but not yet a solution, red square: a fail state (there is no solution), green diamond: an intermediary solution, orange diamond: the optimal solution.)

• all the letters are pairwise dierent: lldifferent(

x S , x E , x N , x D , x M , x O , x T , x Y ) • (1000x S + 100x E + 10x N + x D ) + (1000x M + 100x O + 10x S + x T ) = 10000x M + 1000x O + 100x N + 10x E + x Y • V = 10000x M + 1000x O + 100x N + 10x E + x Y .
The propagation of these constraints leads to a stable state with the domains: , 3, 4, 5, 6, 7, 8} Strategies dene the way to choose variables and for each chosen variable, the way to choose values. If variables are chosen in the order x S , x E , x N , x D , x M , x O , x T , x Y and for each variable, the remaining values in its domain are chosen in an increasing order, the search tree has 29 states with 7 intermediate solutions (a solution which is better than the precedent one). However, if variables are chosen in the order x S , x T , x Y , x N , x D , x E , x M , x O , the search tree has only 13 states with 2 intermediate solutions. These two search trees are shown in Figure 2.6. For each stable state, the left branch corresponds to the case where the chosen variable is assigned to the chosen value, the right branch is the other case where the chosen value is removed from the domain of the chosen variable. In these search trees, blue circle is a stable state but not yet a solution, red square is a fail state (there is no solution), green diamond is an intermediary solution and the orange diamond is the optimal solution. It is worth to notice that it is an exhaustive search, the returned solution is guaranteed to be optimal. Let us notice that in problems where there are several solutions with the same optimal value, the system returns only a single optimal solution.

Dom(x S ) = {9}, Dom(x E ) = {2, 3, 4, 5, 6, 7}, Dom(x M ) = {1}, Dom(x O ) = {0}, Dom(x N ) = {3, 4, 5, 6, 7, 8}, Dom(x D ) = Dom(x T ) = Dom(x Y ) = {2

20

Chapter 2. Constraint Programming

Symmetries in CP

Denition 18 For any CSP instance P = (X, D, C), a solution symmetry σ of P is a permutation of the set X × D that preserves the set of solutions to P . [START_REF] Cohen | [END_REF] In other words, a solution symmetry is a bijection mapping on the set of variable-value pairs, that maps solutions to solutions and also non-solutions to non-solutions. As a result, given a solution symmetry σ, a complete assignment A is a solution if and only if σ(A) is a solution. The solution symmetry makes some sub-trees in the search space to be equivalent. If two sub-trees are equivalent, they will both have no solution, or every solution in a sub-tree corresponds with a symmetric counterpart in the other sub-tree. As a result, it is not necessary to explore both sub-trees which are symmetrically equivalent. In practice, the most common symmetry types are: variable symmetry and value symmetry. Variable (value) symmetries are symmetries that can be represented by a permutation of only variables (values).

Example 16 The 8-queens problem: Find how to place 8 queens on a 8×8 board, so that they do not attack each other: no two queens share the same row, column or diagonal.

This problem can be modelled as a CSP with 8 variables x 1 . . . x 8 for queens, with the domain: dom(x i ) = {1, . . . , 8}. The variable x i represents the column position for the queen in row i. With those denitions of variables, no two queens can be on the same row. The constraints for the problem are:

• ∀i < j ∈ {1, . . . , n} : x i = x j : Those constraints guarantee that two queens cannot be on the same column.

• ∀i < j ∈ {1, . . . , n} : |x i -x j | = j -i: Those constraints guarantee that no two queens are on the same diagonal.

With this modelling, there are solution symmetries as following:

1. The variable symmetry (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 , x 8 ) ↔ (x 8 , x 7 , x 6 , x 5 , x 4 , x 3 , x 2 , x 1 ) that maps x 1 to x 8 , x 2 to x 7 , . . .. It represents a ipping horizontal on the board, as illustrated by the solution 0 and 1 in Figure 2.5. The solution 0 is:

x 1 = 1, x 2 = 5, x 3 , = 8, x 4 = 6, x 5 = 3, x 6 = 7, x 7 = 2, x 8 = 4.
By switching values of x 1 and x 8 , x 2 and x 7 ,. . . , we have the solution 1, which is symmetry to the solution 0:

x 1 = 4, x 2 = 2, x 3 , = 7, x 4 = 3, x 5 = 6, x 6 = 8, x 7 = 5, x 8 = 1.
2. The value symmetry (1, 2, 3, 4, 5, 6, 7, 8) ↔ (8, 7, 6, 5, 4, 3, 2, 1) that maps value 1 to 8, 2 to 7,. . . It represents a ipping vertical on the board, as illustrated by the solutions 0 and 2 in Figure 2.5.

3. The variable-value symmetry that maps every x i = j to x j = i that represents a ipping diagonal on the board, for example the solutions 0 and 3 in Figure 2.5. Another variable-value symmetry maps every x i = j to x 8-j+1 = 8 -i + 1, for instance the solution 0 and 4 in Figure 2.5. Many researches have addressed the importance of symmetry breaking in the search of CSPs. Solution symmetries can be avoided either by adding new constraints to prohibit them or by modifying the search heuristic to prune symmetric states during the search.

Summary

In this chapter, we present basic preliminaries of Constraint Programming. Constraint Programming is a powerful paradigm in which users specify declaratively the problem by variables and constraint, then the solver nds the solutions by constraint propagation and search.

More details on Constraint Programming can be found in text books [Apt 2003] and [START_REF] Rossi | [END_REF] Homogeneity means that the objects inside the same cluster must be similar while well separation expresses that objects in a cluster must be dierent from the objects belonging to other clusters. The problem has been studied since 18 th century with applications in natural sciences, economics, sociology, geology and many other elds. The cluster analysis literature is vast and heterogeneous with hundreds of dierent algorithms.

Clustering could be seen as an optimization problem, i.e., nding a partition of the objects that optimizes a given criterion. In this chapter, we present the state of the art in clustering, mainly from an optimization viewpoint. We focus more on heuristics and exact methods with optimization criteria that are closer to our research. We also bring an overview of constrained clustering with several well-known approaches.

Problem denitions

We consider a dataset of n objects O = {o 1 , . . . , o n }. Each object is described by its values on p attributes, also called variables. We denote the value of the j-th attribute of an object o i by o ij . Most cluster algorithms do not rely on the description of the objects but on a dissimilarity between them.

Suppose that we have a dissimilarity measure between any two objects o i , o j ∈ O, denoted by d ij . The dissimilarity d ij is typically calculated by using a distance metric dened on the attribute space. A distance metric must satisfy the following properties [Han et al. 2011]:

1. Non-negativity:

d ij ≥ 0, ∀i, j ∈ [1, n]. 1 2. Identity of indiscernible: d ii = 0, ∀i ∈ [1, n].

Symmetry:

d ij = d ji , ∀i, j ∈ [1, n].

Triangle inequality:

d ij ≤ d it + d tj , ∀i, j, t ∈ [1, n].
Note that the dissimilarity measure must satisfy the rst three properties, but may not satisfy the triangle inequality. In cluster analysis, the most popular distance metric is Euclidean distance and the squared Euclidean distance [Han et al. 2011]. Given two object o i , o j , the Euclidean distance between them is dened as:

d ij = (o i1 -o j1 ) 2 + (o i2 -o j2 ) 2 + . . . + (o ip -o jp ) 2
The squared Euclidean distance which is dened as:

d ij = (o i1 -o j1 ) 2 + (o i2 -o j2 ) 2 + . . . + (o ip -o jp ) 2
Another popular distance measure is the Manhattan distance [Han et al. 2011], dened as:

d ij = |o i1 -o j1 | + |o i2 -o j2 | + . . . + |o ip -o jp |
Similarly, we denote a similarity measure between two objects o i and o j by s ij . The similarity measure is widely used in spectral clustering. It is typically calculated by a Gaussian function of the dissimilarity [Luxburg 2007]:

s ij = exp(- d ij 2σ 2 )
where σ is a parameter. There exist other similarity measures, for instance the normalized Pearson correlation, the Jaccard measure, the dice coecient measure,. . . [Han et al. 2011].

Most clustering methods can be classied into two categories: hierarchical methods and partitioning methods. A partitioning algorithm aims at nding a partition P while a hierarchical clustering aims at nding a set of partitions P 1 , P 2 , . . . , P l of O.

Criteria 25 Denition 19 Let

C 1 , C 2 , . . . , C k be subsets of O. {C 1 , C 2 , . . . , C k } is a partition of O into k clusters if: • for all c ∈ {1, . . . , k}, C c = ∅ • ∪ c C c = O • for all c = c , C c ∩ C c = ∅.
Many partitioning-based clustering methods rely on the choice of a representative for each cluster, it could be a centroid or a medoid. When the objects are vectors of p quantitative attributes, the centroid of a cluster is computed by averaging the data vectors belonging to that cluster.

Denition 20 The centroid m c of a cluster C c is a point dened as:

∀j ∈ {1, . . . , p} (m c ) j = 1 |C c | o i ∈Cc o ij
When attributes of objects are qualitative, the centroid point can no longer be calculated and a point, usually a medoid, is chosen to represent each cluster.

Denition 21 A medoid x c of a cluster C c is an object of the cluster whose average dissimilarity to all the objects in the cluster is minimal:

x c = arg min o i ∈Cc 1 |C c | x j ∈Cc d ij
Let us notice that a centroid may not belong to O whereas a medoid belongs to O.

Criteria

Criteria for a cluster

Cluster analysis aims at nding clusters, which are homogeneous and well-separated. These requirements are usually expressed by an optimization criterion and the clustering task can be dened as nding a partition of objects that optimizes a given criterion.

The separation of a cluster C c is typically measured by the split criterion.

• the split of a cluster split(C c ) is the minimum dissimilarity between an object of cluster C c and objects belonging to other clusters [Hansen & Jaumard 1997]:

split(C c ) = min o i ∈Cc,o j / ∈Cc d ij
Another criterion for the separation is the cut between clusters. This criterion is widely used in the spectral clustering methods and is measured with the similarities between objects.

• the cut of a cluster C c is the sum of similarities between objects of this cluster and objects belonging to other clusters:

cut(C c ) = o i ∈Cc o j / ∈Cc s ij
The homogeneity of a cluster C c can be measured by the following criteria [Hansen & Jaumard 1997]:

• the diameter diameter(C c ) of cluster C c is the maximum dissimilarity between two of its objects:

diameter(C c ) = max o i ,o j ∈Cc d ij
• the radius radius(C c ) of a cluster C c is measured by the minimum for all objects o i in the clusters of the maximum dissimilarity between o i and any other object

o j in C c : radius(C c ) = min o i ∈Cc max o j ∈Cc d ij
• the within cluster sum of dissimilarities wcsd(C c ) of a cluster C c is measured by the sum of dissimilarities between any two of its objects:

wcsd(C c ) = o i ,o j ∈Cc d ij
Let us notice that, for this criterion the dissimilarities are usually measured by the squared Euclidean distance.

If the objects o i are in the p-dimensional Euclidean space, the homogeneity of a cluster C c can also be measured by:

• the within cluster sum of squares wcss(C c ) of a cluster C c is the sum of the squared Euclidean distances between each object with the centroid m c of the cluster:

wcss(C c ) = o i ∈Cc ||o i -m c || 2
Let us notice that, when squared Euclidean distance is used for measuring the dissimilarities, the wcsd criterion once standardized via the division by the size of each group, is mathematically equivalent to the wcss:

wcss(C c ) = wcsd(C c ) |C c |
• the variance variance(C c ) of a cluster C c is the within cluster sum of squares divided by the cardinality of C c :

variance(C c ) = wcss(C c ) |C c |

Criteria for a partition

The separation of a partition P can be measured either by the minimum separation of all clusters C c with c ∈ [1, n] or by the sum of the separation. This leads to these criteria:

• the split of a partition:

split(P ) = min c∈[1,k] split(C c )
• the sum of splits of a partition:

sum_split(P ) = k c=1 split(C c )
Another criterion that is often used in spectral clustering is the min cut, which minimizes the cut between clusters.

• the cut of a partition is the sum of similarities of edges between cluster:

cut(P ) = 1 2 k c=1 cut(C c )
In case k = 2, the problem of nding min cut is polynomial and can be solved eciently [Stoer & Wagner 1997], otherwise it is NP-Hard. In many cases, the criterion of min cut leads to a partition with imbalanced clusters, for example, one cluster with one vertex and one cluster with the rest. The criterion of ratio cut [Hagen & Kahng 1992] and normalized cut [START_REF] Shi | [END_REF] resolve this problem by introducing the size and the volume of the clusters in the formula.

• the ratio cut of a partition is measured as:

ratio_cut(P ) = 1 2 k c=1 cut(C c ) |C c |
• the normalized cut is measured as:

normalized_cut(C 1 , . . . , C k ) = 1 2 k c=1 cut(C c ) vol(C c ) where vol(C c ) is the sum of similarities in C c vol(C c ) = o i ,o j ∈Cc s ij
The criteria of ratio cut and normalized cut are NP-Hard, even with k = 2.

Chapter 3. Clustering

The homogeneity of a partition P can be measured either by the maximum homogeneity or the sum of homogeneity of each cluster. This leads to the following criteria:

• the maximum diameter:

diameter(P ) = max c∈[1,k] diameter(C c )
• the sum of diameters:

sum_diameter(P ) = k c=1 diameter(C c )
• the within cluster sum of dissimilarities:

wcsd(P ) = k c=1 wcsd(C c )
• the within cluster sum of squares:

wcss(P ) = k c=1 wcss(C c )
Let us notice that, for some criteria, for instance the diameter criterion or the split criterion, this may exist dierent partitions that share the same optimal value. All these criteria lead to NP-Hard problems except the split criterion.

Partitioning methods

Given a dataset of n objects, a partitioning algorithm searches a partition with exactly k (k ≤ n) clusters in which each object belongs to exactly one cluster. In general, a partitioning method optimizes a partitioning criterion. From the point of view of optimization, partitioning methods can be classied into heuristics and exact approaches. In this section, we present well-known partitioning methods that are related to our research.

Heuristics

The k-means and k-medoid algorithms

The k-means algorithm [MacQueen 1967] is the most well-known method for clustering for it is simple and eective. The algorithm is initialized by randomly picking k objects from O as the k initial cluster centroids m 1 , m 2 ,. . . , m k . Then each object is assigned to the closest centroid, forming a partition of the dataset. After that, the cluster centroids are updated and objects are reassigned again to the closest centroid. This process iterates until there is no change of assignments of objects, or equivalently, until there is no change of centroids. The k-means algorithm always converges to a solution where the centroids are stable.

P = {C 1 , C 2 , . . . , C k }
The k-means is ecient if clusters are well separated from each other. It can handle large datasets because the complexity of the algorithm is O(nkt) where t is the number of iterations. However, this method can be applied only when the mean value of clusters can be computed. Moreover, this method is sensible to noises and outliers because those values aect the mean value of clusters. To reduce those eects, the k-medoid algorithm is proposed in [Kaufman & Rousseeuw 1987]. In this method, an actual object is chosen in each cluster as a representative object. This diers from k-means where the representative of a cluster is computed as the mean value of the objects of the cluster. It starts by choosing k objects as the initial cluster representatives. Each object is assigned to its nearest representative object. It then chooses a new representative object for each cluster and the method reassigned the objects by taking into account the new representative. The algorithm iterates until each representative object is actually the medoid. It is more robust to noise and outliers as compared to k-means.

The Furthest Point First algorithm

The method FPF (Furthest Point First) introduced in [Gonzalez 1985] relies on the maximum diameter criterion. In this method, each cluster has an object as the head and every object in a cluster satises the property of being closer to their head than to the other heads. The method starts by choosing an object as the head of the rst cluster and assigns all objects to it. In the next iteration, the furthest object from the rst head is chosen as the head of the second cluster. Any object which is closer to the second head than to the rst one is moved to the second cluster. The algorithm iterates: choosing the object that is the furthest to its head as the head of the new cluster and reassigning objects. It stops after k iterations since k clusters are needed. The time complexity is O(kn). Let d opt be the global optimal, then this heuristic is guaranteed to nd a clustering with a maximum diameter d such that d opt ≤ d ≤ 2d opt if the dissimilarity measure satises the triangle inequality. Moreover, [Gonzalez 1985] showed that this approximation ratio is tight: nding d such that d ≤ (2 -)d opt is NP-hard for all > 0 when objects are in a 1 Randomly select an object, mark it as head 1 , the head of cluster 1

2 foreach i from 2 to k do 3
Select the object that is the furthest to its head, mark it as head i 4

Assign every object o to cluster i if o is closer to head i than to its head 5 end

Spectral methods

The main idea of spectral methods is to change the presentation of n objects x i to y i by the Laplacian graph. Then the methods use the classic k-means clustering algorithm to partition new objects y i . Spectral methods use the graph of similarity G = (V, E) where each vertex v i represents an object o i , two vertices are connected if the similarity between them is below a given threshold. G is a weighted graph where the weight of edges v i and v j , denoted w ij , is the similarity s ij between o i and o j . The graph G is undirected and the degree of a vertex v i is measured as the sum of weights of the edges connected to v i :

degree(i) = n j=1 w ij
The degree matrix D is a diagonal matrix, where the diagonal contains the degrees of each

vertex v i of G D ii = degree(i), D ij = 0 ∀i = j
The spectral methods are used to solve the clustering problem with the criteria of min cut, ratio cut or normalized cut. The optimization problems with those criteria are NP-Hard. However, they can be relaxed by using spectral concepts of graph analysis. The relaxation is done by the Laplacian graph.

The unnormalized Laplacian graph is introduced in [Chung 1997] :

L = D -W
where W is the matrix of (w ij ).

The unnormalized Laplacian graph has the following properties:

• For every vector f = (f 1 , . . . , f n ) ∈ R n : f t Lf = 1 2 n i=1 n j=1 w ij (f i -f j ) 2
• L is symmetric and positive semi-denite

• The smallest eigenvalue of L is 0, it corresponds to the constant one vector.

• L has n non-negative eigenvalues

0 = λ 1 ≤ λ 2 ≤ . . . ≤ λ n
In the literature, there are dierent types of Laplacian graph and each Laplacian graph has special properties, suitable for dierent criteria. The basic spectral clustering algorithm is presented in Algorithm 5.

Algorithm 5: Spectral methods 1 Compute the similarities between objects s ij , ∀i, j ∈ [1, n] 2 Compute the Laplacian graph 3 Compute the rst k eigenvectors u 1 , . . . , u k of the Laplacian graph 4 Construct a matrix U ∈ R n×k by those k eigenvectors 5 The new dataset consists of n objects, y i is the vector corresponding to the i-th row of U 6 Apply the k-means method for the new n objects y i [Luxburg 2007] shows how the spectral method works with the criterion of minimizing the ratio cut, in the case k = 2, with the unnormalized Laplacian graph .

For a cluster C c , we dene the vector f = (f 1 , f 2 , . . . , f n ) t with:

f i = |O \ C c | |C c | , ∀i : o i ∈ C c and f j = - |C c | |O \ C c | , ∀j : o j / ∈ C c (3.1)
Due to the rst property of the Laplacian graph , we have:

f t Lf = 1 2 n i=1 n j=1 w ij (f i -f j ) 2 = 1 2 o i ∈Cc o j / ∈Cc w ij |O \ C c | |C c | + |C c | |O \ C c | 2 + 1 2 o i / ∈Cc o j ∈Cc w ij - |O \ C c | |C c | - |C c | |O \ C c | 2 = o i ∈Cc o j / ∈Cc w ij |O \ C c | |C c | + |C c | |O \ C c | 2 = cut(C c ) |O \ C c | |C c | + |C c | |O \ C c | + 2 = cut(C c ) |O \ C c | + |C c | |C c | + |O \ C c | + |C c | |O \ C c | = n × ratio_cut(C c )
We also have:

n i=1 f i = o i ∈Cc |O \ C c | |C c | - o i / ∈Cc |C c | |O \ C c | = |C c | |O \ C c | |C c | -|O \ C c | |C c | |O \ C c | = 0
The vector f is orthogonal to the constant vector 1, and it satises:

||f || 2 = n i=1 f 2 i = |C c | |O \ C c | |C c | + |O \ C c | |C c | |O \ C c | = n
The problem of minimizing the ratio cut is equivalent with:

min Cc f t Lf subject to: f dened in equation 3.1, f ⊥1, ||f || = √ n
This is a discrete optimization problem and it is NP-hard. This problem can be relaxed allowing f i to take arbitrary values in R. The problem becomes a relaxed optimization problem:

min Cc f t Lf subject to: f ⊥1, ||f || = √ n
The solution of this relaxed optimization problem is given by the eigenvector of the second smallest eigenvalue of L. To obtain the partition corresponding to the approximate ratio cut, the vector f must be transformed into a discrete vector. A simple way for the transformation is:

o i ∈ C c if f i ≥ 0, o i / ∈ C c if f i < 0
In the case k > 2, a better way for the transformation is to consider f i as points and then to cluster them by k-means algorithm (as presented in Algorithm 5).

The main advantages of spectral methods are that they are not aected by particular forms of clusters whereas other methods, for example, k-means algorithm tends to nd convex sets. It nds clusters based on the connectivity, not on the compactness like k-means. It can solve rather large datasets because it solves a linear problem.

Exact approaches

Graph-theoretical algorithm

Graphs are formed by a set of vertices and a set of edges, connecting pairs of vertices. If we consider objects as vertices and using the dissimilarities as the edges between them, the clustering analysis is close to the graph theory. The problem of nding an optimal partition with the criterion of maximum diameter is proved to be NP-Hard and an exact method based on the problem of graph colouring is proposed in [START_REF] Hansen | [END_REF], presented in Algorithm 6. The algorithm is based on the observation that: there are at most n(n -1)/2 distinct values of dissimilarities. Suppose that we have a partition P with k clusters and d max is the maximum diameter. It means that every pair of objects (o i , o j ) in which d ij > d max must be in dierent groups. Let G = (V, E) be a graph where V has n vertices corresponding to the n objects o 1 , o 2 , . . . , o n , and E contains edges

(o i , o j ) satisfying d ij > d max .
It is evident that the graph G can be coloured into k colours. Indeed, from the partition P, we can use the same colour for all objects of the same cluster, thus we need at most k colours for k clusters.

Algorithm 6: graph colouring approach for clustering with the diameter criterion

1 Compute the matrix of dissimilarity D = (d ij ) 2 Build the initial graph G = (V, E) where V contains n objects and E is empty 3 d max ← max i,j∈[1,n] d ij 4 repeat 5 Add edges (o i , o j ) in which d ij = d max to E

Branch-and-bound algorithms

As the problem is dened as assigning n objects into k dierent clusters, with a criterion to optimize, branch-and-bound approaches can be useful to solve the problem. The performance of a branch-and-bound method depends on how the algorithm calculates the bounds and on the branching strategy. A branch-and-bound algorithm for the criterion of the within cluster sum of squares has been proposed in [Koontz et al. 1975]. It was improved later in [Brusco 2006]. For the criterion of within cluster sum of dissimilarities, a branch-and-bound algorithm has been proposed in [Klein & Aronson 1991], and in [Brusco & Stahl 2005]. For the criterion of minimum diameter, a method has been proposed in [Brusco 2003].

The method of branch-and-bound is based on the evaluation of partial (incomplete) partitions. A partial partition Q is an assignment of p objects to k clusters and it remains n -p unassigned objects. It may exist empty clusters in a partial partition.

Depending on the criteria, dierent bounds can be used to evaluate the partial solution. The basic BAB algorithm which uses a depth-rst-search strategy is presented in Algorithm 7, supposing that the objective is to minimize a criterion.

Return

The real challenge is how to eciently evaluate the partial solutions in order to prune useless branches in the search tree. To determine whether a partial solution may be a part of the optimal solution, a typical way is to compute a lower bound of the criterion. If the lower bound is greater than the upper bound, then the partial solution cannot be a part of the optimal solution and the search backtracks to explore other branches.

For the criterion of diameter [Brusco 2003] uses three tests to evaluate the partial solution:

• The rst test is used to check if there are enough unassigned objects with respect to the remaining empty clusters.

• The second test checks if the diameter of the partial solution is equal to or greater than the maximal diameter of the current best solution, based on assigned objects.

• The third test is based on unassigned objects, it tests if there exists an object that cannot be assigned to a cluster without leading to a new diameter equal to or greater than the maximal diameter of the current best solution.

In any of these three tests return true, there is no need to continue with this partial solution and the search backtracks.

For the criterion of the within cluster sum of dissimilarities The branch-andbound algorithm in [Klein & Aronson 1991] computes a lower bound of wcsd for the partial solutions. The lower bound of wcsd is computed as a sum of three parts: the rst part is the sum of dissimilarities of assigned objects, the second part is the sum of dissimilarities between assigned objects and unassigned ones and the third part is the sum of dissimilarities between unassigned objects.

As the rst part is the sum of dissimilarities of assigned objects, it can be computed exactly. The second and the third part cannot be calculated exactly and a lower bound for these parts are measured by heuristics. [Brusco & Stahl 2005] proposed to improve the lower bound of wcsd by computing the optimal wcsd of unassigned objects for the third part instead of using a heuristic. [Brusco & Stahl 2005] solves the problem in a Repetitive Branch-and-Bound Algorithm (RBBA). Given a dataset of n point o 1 , . . . , o n and a clustering task dividing the points into k clusters, the RBBA is realized in n -k steps:

• Step 1: Find the optimal wcsd of k + 1 objects {o n-k , o n k +1 , . . . , o n } into k clusters
with a branch-and-bound algorithm. The search strategy assigns objects into clusters following the order: o n-k , . . . , o n .

•

Step 2: Find the optimal wcsd of k + 2 objects {o n-k-1 , o n k , . . . , o n }.

• ...

•

Step n -k: Find the optimal wcsd with n objects {o 1 , . . . , o n }.

Let us consider the step i that nds the optimal wcsd of k For the criterion of within cluster sum of squares The RBBA algorithm can be applied to the wcss criterion [Brusco 2006, Carbonneau et al. 2012]. The basic idea is the same. For this criterion, a lower bound of a partial solution can be measured by the sum of wcss of assigned objects and the optimal wcss of unassigned objects. The RBBA solves the problem rst with k + 1 objects, then with k + 2 objects and so on, until all n objects are considered.

+ i objects {o n-k-i+1 , . . . , o n } into k

Hierarchical clustering methods

Hierarchical clustering aims at nding a hierarchy of clusters. Unlike partitioning methods, hierarchical clustering nds a sequence of nested partitions. The result is a tree diagram which is called a dendrogram. There are two types of hierarchical clustering:

• Agglomerative: a "bottom up" approach, it begins with n clusters, each cluster contains one and only one object. At each iteration, two "close" clusters are merged until a stopping criterion is satised.

• Divisive: a "top down" approach, it begins with only one cluster that contains all the n objects. It successively splits one cluster into two smaller parts until a stopping criterion is satised.

The agglomerative hierarchical clustering algorithms start by a partition P n with n clusters C 1 to C n , each cluster C i contains exactly one object o i . The basic algorithm is presented in Algorithm 8. At each iteration, the algorithm selects two close clusters based on a criterion and merges them until there is only one cluster. The most popular agglomerative hierarchical clustering algorithms are single-linkage and complete-linkage algorithms. The main dierence between those two algorithms is how they calculate the distance between two clusters. In the single-linkage algorithm, the distance between two clusters is the minimum of the distances between any pairs of objects of the two clusters whereas in the complete-linkage algorithm, it is the maximum of the distances of all pairs of objects in the two clusters. The single-link algorithm may suer from the chain eect [Johnson 1967], i.e. a chain of close objects may group very dierent objects in the same cluster whereas the complete-link cluster may produce compact clusters but not well-separated.

It is proved in [Delattre & Hansen 1980] that at all levels, the single-linkage algorithm provides an optimal partition with respect to the criterion of split, so it can be used to nd the exact split with any number of clusters k.

Algorithm 8: Basic agglomerative hierarchical clustering algorithm 1 Construct a partition with n clusters

P n = {C 1 , C 2 , . . . , C n }, C i = {o i } 2 foreach t ∈ [1, n -1] do 3 Select C i , C j ∈ P n-t+1 based on a criterion 4 C n+t = C i ∪ C j 5 P n-t = (P n-t+1 ∪ C n+t ) \ {C i , C j } 6 end
In contrast to the agglomerative algorithms, the basic divisive hierarchical clustering, presented in Algorithm 9 starts by a partition P 1 which contains only one cluster C 1 with the whole objects and it divides a cluster based on a criterion, until there are n clusters.

Algorithm 9: Basic divisive hierarchical clustering algorithm 1 Construct a partition with only one cluster

P 1 = {C 1 }, C 1 = {o 1 , o 2 , . . . , o n } 2 foreach t ∈ [1, n -1] do 3 Select C i ∈ P t based on a criterion 4 Divide C i into C 2t and C 2t+1 5 P t+1 = (P t ∪ C 2t ∪ C 2t+1 ) \ C t 6 end

Constrained Clustering

In order to better model the task, but also in the hope of reducing the complexity, userspecied constraints are added, leading to Constrained Clustering that aims at nding clusters that satisfy user-specied constraints. User constraints can be classied into cluster-level constraints, specifying requirements on clusters, or instance-level constraints, specifying requirements on pairs of objects.

Instance-level constraints (also called pairwise constraints) are constraints on pairs of objects. There are two types of instance-level constraints, introduced for the rst time in [START_REF] Wagsta | [END_REF], the must-link and cannot-link constraints. A must-link constraint between two object o i , o j , denoted by M L(o i , o j ), expresses that both objects o i and o j must be in the same cluster. In contrast, a cannot-link constraint between o i and o j , denoted by CL(o i , o j ), means that these two objects must not be in the same cluster. According to [START_REF] Davidson | [END_REF], the instance-level constraints have the following properties:

• Must-link constraints are transitive: Let CC a and CC b be connected components (completely connected subgraphs by must-link constraints), let o i and o j be objects [START_REF] Davidson | [END_REF] in CC a and CC b respectively, then:

M L(o i , o j ), o i ∈ CC a , o j ∈ CC b ⇒ M L(o x , o y ), ∀o x , o y : o x ∈ CC a , o y ∈ CC b
• Cannot-link constraints can be entailed. Let CC a and CC b be connected components (completely connected subgraphs by must-link constraints), let o i and o j be objects in CC a and CC b respectively, then:

CL(o i , o j ), o i ∈ CC a , o j ∈ CC b ⇒ CL(o x , o y ), ∀o x , o y : o x ∈ CC a , o y ∈ CC b
The idea of must-link and cannot-link constraints might seem simple but in fact, those constraints are powerful in many applications. Many researchers have reported that increasing the number of pairwise constraints may improve the accuracy of the result. The must-link and cannot-link constraints can also be used to express other user-constraints. [Davidson & Ravi 2005b] has introduced -constraint and δ-constraint which can be handled with pairwise constraints:

• δ-constraint (also called minimum split constraint) expresses that the distance between any pair of points which are in two dierent clusters must be at least δ. It can be handled by a conjunction of must-link constraints between all pair of objects with distance less than δ.

• -constraint expresses that: for any cluster C c containing more than one object, for any object o i ∈ C c , there must be another object o j ∈ C c such that the distance between o i and o j does not exceed : d ij ≤ . This constraint tries to capture the notion of density, introduced in DBSCAN [Ester et al. 1996]. It is shown in [Davidson & Ravi 2005b] that this constraint is equivalent with a disjunction of must-link constraints. For each point o i , a set X of point o j such that d ij ≤ is computed, then a disjunction of must-link constraints between o i and points in X are generated.

[ [START_REF] Davidson | [END_REF] show that nding feasibility partition with those constraints is not an easy task. Table 3.1 expresses the complexity of nding a feasibility partition with dierent constraints.

Another constraint that can be expressed by instance-level constraints is the maximum diameter constraint. This constraint species an upper bound γ on the diameter of the clusters to express that any pair of points which are in the same cluster must be at a distance of at most γ. This constraint can be handled by conjunction of cannot-link constraints between all pair of objects with distance greater than γ.

Clustering algorithms for instance-level constraints

In the last ten years, many works have been done to extend classical algorithms for handling must-link and cannot-link constraints, as for instance an extension of COBWEB [START_REF] Wagsta | [END_REF], k-means [START_REF] Wagsta | [END_REF], Bilenko et al. 2004], hierarchical non supervised clustering [Davidson & Ravi 2005a] or spectral clustering [START_REF] Lu & Carreira-Perpinan | [END_REF][START_REF] Wang | [END_REF], . . . Those methods can be classied into two general approaches [START_REF] Davidson | [END_REF]: constraint-based and distance-based methods. In constraint-based approaches, the clustering algorithm is modied to integrate pairwise constraints whereas in distance-based approaches, only the distance measure is modied.

3.5.1.1 Constraint-based methods [START_REF] Wagsta | [END_REF] proposed a modied k-means algorithm, called COP-Kmeans to integrate must-link and cannot-link constraints. The major modication is that, at each iteration, objects are updated so that none of the constraints is violated. The drawback of this approach is that it tries to satisfy all the constraints but the algorithm does not provide any backtracking technique. As a result, the algorithm may fail to nd a partition even if one exists. It happens when there are many constraints, especially cannot-link constraints.

Algorithm 10: COP-Kmeans algorithm 1 Generate randomly k objects as initial cluster centroids

o j : o j ∈ C c ∩ CL(o i , o j ) o j : o j / ∈ C c ∩ M L(o i , o j ) 5 if no such cluster C c exists then 6 Return FAILED 7 end 8 end 9
Recalculate the new centroid of each cluster as the mean value of the objects in each cluster must-link and cannot-link constraints by modifying the objective function. In K-means algorithm, at each iteration, the centroid m c is updated as:

m c = o i ∈Cc o i |C c |
After that, each object is reassigned to the closest centroid to minimize the vector quantization error (VQE):

V QE = 1 2 c∈[1,k] o i ∈Cc ||m c -o i || 2
The CVQE algorithm modies the error function to penalize the violated constraints, [Davidson & Ravi 2005b] denes the constrained vector quantization error CV QE c of a cluster C c by:

CV QE c = 1 2 o i ∈Cc (m c -o i ) 2 + 1 2 o i ∈Cc,M L(o i ,o j ),o j ∈C c ,c =c ||m c -m c || 2 + 1 2 o i ∈Cc,CL(o i ,o j ),o j ∈Cc ||m c -m h(c) || 2
where h(c) returns index of the cluster (other than c) whose centroid is the closest to m c . In this equation, the rst part is the VQE criterion as for the k-means algorithm. The second part is the penalty for the violation of must-link constraints. If a constraint M L(o i , o j ) is violated, the penalty is measured as the square distance between the cluster centroids containing objects of the constraint. The third part is the penalty for cannot-link constraints, it is measured as the squared distances between the centroid containing both objects and the nearest centroid.

The objects that appear in no constraint are always assigned to the nearest centroid. For the other objects, they are assigned to the centroid so as to minimize the CVQE.

The update rule of centroids are:

m c = o i ∈Cc o i + M L(o i ,o j ),o j ∈C c ,c =c m c + CL(o i ,o j ),o j ∈Cc m h(c) |C c | + M L(o i ,o j ),o j ∈C c ,c =c 1 + CL(o i ,o j ),o j ∈Cc 1
The meaning of the updating on centroids is as follows: if a must-link constraint M L(o i , o j ) is violated, the centroid m c of the cluster containing o i moves closer to the centroid m c of the cluster containing o j . Similarly, the violation of a cannot-link constraint CL(o i , o j ) makes the centroid containing both constrained objects move to the nearest cluster centroid.

[ [START_REF] Pelleg | [END_REF] proposed LCVQE, an improved algorithm of CVQE, by changing the penalty for violating constraints. For the violation of a constraint M L(o i ∈ C c , o j ∈ C c ), suppose that o i is assigned before o j , the penalty is now the squared distance between o j and m c . For the violation of a constraint CL(o i ∈ C c , o j ∈ C c ), the point that is further to m c is chosen. The penalty is the squared distance between this point to the nearest centroid to m c . The value LCVQE of each cluster c is measured between o i and o j as:

LCV QE c = 1 2 o i ∈Cc (m c -o i ) 2 + 1 2 o i ∈Cc,M L(o i ,o j ),o j / ∈C c ,c =c ||m c -m c || 2 + 1 2 o i ∈Cc,CL(o i ,o j ),o j ∈Cc ||m c -m h(c) || 2
By modifying the formulation, the algorithm LCVQE has better computation time.

Distance-based algorithms

Another way to incorporate instance-level constraints is to modify the distance metric or the objective function. The basic idea is that if there is a must link constraint between two objects o i and o j , the distance between them might be smaller than usual, so they have a higher chance to be in the same cluster, similarly to the cannot-link constraints.

[ Xing et al. 2003] proposed an algorithm to learn a distance metric over must-link and cannot-link constraints. The distance metric learning is formulated as an optimization problem. The objective function is the sum of distances of pairs in the must-link constraints. The algorithm searches to minimize this sum, under the condition that the distances of pairs in the cannot-link constraints are greater than a constant. The distance metric is of the form:

d(o i , o j ) = d A (o i , o j ) = ||o i -o j || A = (o i -o j ) T A(o i -o j )
A is semi-denite positive and if A = I, we have the Euclidean distance. The problem is now to nd the matrix A:

min M L(o i ,o j ) d 2 A (o i , o j ) with CL(o i ,o j ) d A (o i , o j ) ≥ 1 A 0
The optimization problem is convex and dierentiable and it is possible to nd a global optimum in polynomial time.

Clustering algorithms for cluster-level constraints

Cluster-level constraints specify requirements on clusters. Beside the -constraint (minimum split constraint), δ-constraint and maximum diameter constraints, other examples of cluster-level constraints are:

• The minimum capacity constraint requires that each cluster has a number of objects greater than a given threshold α:

∀c ∈ [1, k], |C c | ≥ α
• The maximum capacity constraint requires each cluster to have a number of objects inferior to a predened threshold

β: ∀c ∈ [1, k], |C c | ≤ β.
• The balance constraint expresses that all clusters have approximately the same size, or the fraction between the size of smallest and largest cluster must be greater than a given threshold θ: 2006] proposed a scalable framework for clustering to incorporate the balance constraint. The framework consists of three steps:

min i∈[1,k] |C i | max j∈[1,k] |C j | ≥ θ [Banerjee & Ghosh
• sampling the data into a small representative subset of the points,

• clustering of the sampled data,

• populating the initial clusters with the remaining data.

The rst step is for the scalability of the framework. In the second step, any clustering algorithm can be used. The third step consists of populating and rening the clusters. In populating phase, the remaining objects are assigned to the existing clusters such that the constraint is satised while ensuring good quality clusters. In rening phase, objects are re-assigned to improve the objective function.

[ Demiriz et al. 2008] proposed an approach to integrate the minimum capacity constraint to k-means algorithm, to avoid empty clusters or clusters with few points in the result. [Demiriz et al. 2008] modies the cluster assignment step of k-means algorithm to introduce the minimum capacity constraint, then solves the constrained problem by using linear programming or network simplex methods.

Summary

In this chapter, we made a brief introduction to clustering. We presented important algorithms, mainly from an optimization viewpoint, with a focus on those related to our research. Afterwards, we presented basic notions of constrained clustering and several well-known approaches. Most researches aim at highly optimized and scalable algorithms, tailored towards specic tasks, in imperative paradigm. Researchers express algorithms in a way that machine understand how to solve the problems. In return, they expect that the machine will nd out what they want. Recently, there has been a growing interest in declarative programming for data mining [START_REF] Raedt | [END_REF]. According to [START_REF] Lloyd | [END_REF]]: "A declarative programming involves stating what is to be computed, but not necessarily how it is to be computed". Here, researchers tell the machine what they want and the machine will gure out how to do it. Although declarative framework may be less ecient than tailored algorithms for a specic task, a declarative framework is way more exible, easy to incorporate new knowledge to the task [START_REF] Raedt | [END_REF]. Table 4.1 presents a comparison between Declarative and Imperative approach.

To solve a problem in a declarative paradigm, users specify it in a modelling language and a solver automatically generates the solutions. Declarative frameworks can be expressed in Constraint Programming, Propositional Satisability, Integer Linear Programming or Logic Programming. The real challenge is how to encode the problem eectively. In this Chapter, we present state of the art of recent declarative frameworks for Data Mining, specially researches that are related to cluster analysis.

Constraint Programming framework for k-pattern set mining

Pattern mining is a data mining method that involves nding all existing patterns of a database that satisfy a set of constraints. A well-known example of pattern mining is nding frequent itemsets in market basket analysis. We consider each product is an item and a group of items is an itemset. Transactions are instances of products that are bought by each customer. An itemset, such as milk and bread, is frequent if these products appear frequently together in a transaction dataset.

[De [START_REF] Raedt | Constraint programming for itemset mining[END_REF]] introduces a framework in Constraint Programming for itemset mining, which is exible and allows to combine dierent mining constraints.

[ [START_REF] Khiari | [END_REF] proposes a generic approach to express and mine n-ary patterns, i.e. partterns involves several local pattern. [Guns et al. 2013] introduces the problem of kpattern set mining that nds a set of k related patterns under constraints. The set of constraints consists of local constraints, i.e. constraint on individual pattern, and global constraints, i.e. constraints on the overall pattern set. [Guns et al. 2013] proposes a framework in Constraint Programming for k-pattern set mining. The framework is similar to [START_REF] Khiari | [END_REF]], but it covers wider range of tasks.

The framework can be applied in various well-known mining tasks, including conceptual clustering.

The principle of conceptual clustering is not limited to nd a partition of the data, but it has to generate a conceptual description for each cluster. A cluster concept is a set of properties that every object of the cluster has. If objects are transactions, the description of a cluster can be represented by a set of items and the conceptual clustering can be considered as a k-pattern set mining problem under specic constraints.

To solve the problem of k-pattern set mining, [Guns et al. 2013] proposes a general framework in Constraint Programming. Considering a dataset of n items and m transactions, the goal is to nd k-patterns π 1 , . . . , π k that satisfy all local and global constraints. The model uses k × n binary variables I[p, i] ∈ {0, 1} for items and k × m variables T [p, t] ∈ {0, 1} for transactions (p ∈ {1, . . . , k}, i ∈ {1, . . . , n}, t ∈ {1, . . . , m}). The variable I[p, i] = 1 indicates that item i is in the pattern π p . Similarly, the variable T [p, t] = 1 indicates that transaction t is in the pattern π p . By using Constraint Programming, the model is exible to express dierent kinds of constraints on items, on transactions, on patterns and on an objective function:

• Coverage constraint: the k patterns are disjoint.

• The k patterns cover all transactions.

• Each pattern is closed: the itemset must be the largest itemset that is contained in all transactions of the pattern.

• Constraint on the size of the pattern which is measured by the size of the itemset.

• Constraint on the dierence between the sizes of the patterns.

• Constraint on the number of transactions in each pattern.

• Constraint on the overlap: it allows a percentage of transactions that can be in more than one pattern/cluster.

The model integrates dierent objective functions: maximizing the minimum cluster size or minimizing the dierence between cluster sizes. Although the proposed CP model, when applied to cluster analysis, is limited to transaction database and conceptual clustering, the model is exible in modeling constraints and objective functions.

SAT based framework for clustering

In propositional logic, a propositional variable x i is a boolean variable, which takes a truth value 0 (false) or 1 (true). A literal l i is a variable x i or its negation x i , a clause C i is a disjunction of literals l i , for example:

C 1 = x 1 ∨ x 2 ∨ x 3 . A conjunctive normal form (CNF) formula is a conjunction of clauses C i , for example: C 1 ∧ C 2 ∧ C 3 . A literal x i (resp.
x i ) is satised if x i gets value 1 (resp. 0). A clause is satised if at least one of its literals is satised, and a CNF is satised if all the clauses are satised. The SAT problem consists of nding an assignment of truth values to variables that satises a CNF formula. A 2-SAT problem is a SAT problem in which each clause involves no more than 2 variables. In MaxSAT problem, the objective is to satisfy the maximum number of clauses. A partial MaxSAT instance is a CNF formula which contains soft (relaxable) and hard (non-relaxable) clauses. Solving a partial MaxSAT problem consists of nding an assignment that satises all hard clauses and the maximum number of soft clauses. The SAT problem was the rst known example of an NP-complete problem (except for the 2-SAT problem, which has a polynomial time algorithm). However, modern SAT solvers can eciently solve problems containing a large number of clauses.

A Constraint Language based on SAT for Declarative Pattern

Mining [Metivier et al. 2011] presents a constraint-based language that can be used for modelling dierent pattern mining problems. The general idea is similar to the framework of [Guns et al. 2013], but the resolution is completely dierent. A pattern mining problem can be dened as queries of the language. A query is a conjunction of constraints Chapter 4. Declarative Data Mining and all primitive constraints of the language are modelled using the SAT framework. The problem then can be solved by a SAT solver. In [START_REF] Métivier | [END_REF], the authors illustrate how to use the language for the problem of clustering on transaction datasets. Each cluster is a pattern that consists of a set of transactions. The problem of nding k patterns/clusters C 1 , . . . , C k can be expressed by the conjunction of the following constraints in the language:

• ∧ i∈[1,k] isNotEmpty(C i ): each pattern must contain at least one transaction.

• coverTransaction(C 1 , . . . , C k ): each transaction is covered by at least one pattern.

• noOverlapTransactions(C 1 , . . . , C k ): there is no transaction, which is covered by more than one pattern.

• canonical(C 1 , . . . , C k ): the pattern is in lexicographic order. This constraint is for symmetry breaking.

The language also supports must-link, cannot-link and size constraints, which allow users to handle dierent constrained clustering tasks, by changing or adding constraints to the modelling. A SAT solver then can nd some or all partitions that satisfy all the queries. The SAT encoding of the language integrates features of SAT solvers to improve the performance, for example: binary clauses, unit propagation, sorting networks. The experiments show that with this encoding, SAT solvers can solve eciently the clustering task with medium size datasets. For example, with the MiniSat solver, the rst ten solutions are found quickly for the dataset Mushroom with 8124 transactions. In conclusion, this constraint-based language is declarative and exible, but when applied to constrained clustering, it is limited to itemset data and does not integrate optimization criteria.

4.2.2 A SAT framework for constrained clustering with 2 clusters [Davidson et al. 2010] proposed a 2-SAT framework for constrained clustering, allowing nding a partition with 2 clusters that optimizes the diameter criterion or the split criterion. For the modelling, [Davidson et al. 2010] uses n boolean variables x i ; x i is assigned to 1 (0) if the object o i is aected to cluster 1 (0). The assignment of truth values to these n variables gives us a partition that consists of two clusters. With those variables, the instance level constraints can be easily expressed by CNF formulas. A must-link constraint between objects o i and o j is expressed by:

(x i ∨ x j ) ∧ (x i ∨ x j )
Similarly, a cannot-link constraint between objects o i and o j is expressed by :

(x i ∨ x j ) ∧ (x i ∨ x j )
The diameter constraint α can be handled by adding cannot-link constraints to every pair of points o i and o j such that their distance exceeds α. With this modelling, a diameter constraint is expressed by adding two clauses (x i ∨ x j ) ∧ (x i ∨ x j ) for every pair o i and o j such that d ij > α. Similarly, a split constraint β is handled by adding two clauses (x i ∨ x j ) ∧ (x i ∨ x j ) to every pair o i and o j such that d ij ≤ β. With these CNF formula, a SAT solver can be used to nd a partition with 2 clusters that satises the must-link, cannot-link, diameter and split constraints. To nd a partition that optimizes the diameter criterion (resp. the split criterion), [Davidson et al. 2010] uses repeated calls to solve the SAT formula. It is clear that the optimal value of diameter (split) criterion must be one of the pairwise dissimilarities. Therefore, the optimal value can be found by a binary search. For example, with the diameter criterion, at each iteration, a dissimilarity x is chosen for a constraint on the diameter (diameter ≤ x), then a SAT formula is generated. Then a SAT solver is used to nd the assignment of truth values to variables. The assignment gives us a partition that satises all the constraints, including the diameter constraint. If such a solution exists, the optimal diameter is less or equal to the dissimilarity x, otherwise the optimal diameter is greater than x. Suppose that there are m dierent dissimilarities, the algorithm takes at most O(log m ) calls to nd the optimal diameter.

The approach is somewhat similar to the algorithm of [START_REF] Hansen | [END_REF]]:

[ Davidson et al. 2010] uses a SAT formula to test if a dissimilarity can be the optimal diameter while [START_REF] Hansen | [END_REF]] uses a graph coloring method to verify that. This framework is the rst modelling of a clustering problem as an instance of SAT. It is limited to 2 clusters and the optimization criterion is not modelled but is expressed by the search strategy, which repeatedly generates the SAT formula and then calls a SAT solver.

MaxSAT framework for Optimal Constrained Correlation Clustering

Correlation clustering [Bansal et al. 2004] is a clustering method that nds the partition based on the relationships between the objects, instead of the actual representations of the objects. Correlation clustering uses a boolean similarity measure s ij between two objects o i and o j : s ij = 1 if o i and o j are similar and s ij = 0 if they are dissimilar. In correlation clustering, the objective is to minimize the cost function:

i,j:λ i =λ j (1 -s ij ) + i,j:λ i =λ j s ij
where λ i denotes the index of the cluster the object o i is assigned to. The goal of correlation clustering is to nd the partition that agrees as much as possible with the similarity measure. This criterion does not require to specify the number of clusters k as a parameter, as the optimal number of clusters could be any value between 1 and n [Bansal et al. 2004]. Finding a partition that optimizes this criterion is NP-hard. [START_REF] Berg | [END_REF] presents a partial MaxSAT framework for nding a global optimal solution for Correlation Clustering. The Correlation Clustering is encoded as a partial MaxSAT instance, which contains a set of hard clauses and a set of soft clauses. With this encoding, a MaxSAT solver is used to nd an optimal solution of the instance. This solution can be mapped to the optimal solution of the Correlation Clustering. Instancelevel constraints can also be encoded in the MaxSAT instance.

The encoding uses boolean variables x ij , where i < j, i, j ∈ {1, . . . , n}. x ij is equal to 1 if the two objects o i and o j are in the same cluster. Hard clauses are used to ensure a well-dened clustering: for each triple of objects (o i , o j , o t ), if o i , o j are in the same cluster and o j , o t are in the same cluster, then objects o i , o t must be in the same cluster. The condition is equivalent with: if x ij = 1 and x jt = 1 then x it = 1. It can be expressed by a hard clause:

(x ij ∨ x jt ∨ x it )
The soft clauses are used to encode the cost function. To minimize the cost function, each pair of objects (o i , o j ) with the similarity s ij = 1 should be in the same cluster, hence x ij should be equal to 1, or it can be expressed by a soft clause (x ij ). Similarly, with each pair of objects (o i , o j ) that s ij = 1, a soft clause (x ij ) is added to the MaxSAT formula.

With this encoding, the instance-level constraints can be easily expressed by hard clauses. A must-link constraint between objects o i and o j is expressed by a clause (x ij ) and a cannot-link constraint between objects o i and o j is expressed by a clause (x ij ). From a MaxSAT solution, it is easy to construct a clustering as follows:

• Assign object o 1 to the rst cluster C 1 ,

• For every i such that variable x 1i = 1, assign object o i to cluster C 1 ,

• Find the smallest i for which o i is not assigned to a cluster, then assign object o i and any object j such that x ij = 1 to cluster C 2 ,

• Iterate this process until every object is assigned.

It is demonstrated in the paper that the partition found by this method is an optimal partition. In this encoding, the number of clusters k is not bounded and the optimal partition may have a large number of clusters. The authors also propose an alternative MaxSAT formulation, which uses an upper bound k on the number of clusters.

Column generation framework for Constrained Sum of Squares Clustering

Column generation is an ecient method for solving linear programming (LP) problems with a large number of variables. In many LP problems, most of the variables will be zero in the optimal solution and only a subset of variables is needed to solve the problem. In a column generation method, the LP problem is split into the master and the auxiliary problem. The master problem is the initial problem. The column generation method begins by solving a restricted master problem, which is the master problem but with only a small number of variables. The auxiliary problem is a problem created from the solution of the restricted master problem. It is used to nd a new variable with negative reduced cost, i.e. a variable that can improve the objective function of the master problem. If such a variable exists, it is added to be considered for resolving the restricted master problem.

In the column generation method, the restricted master and the auxiliary problems are solved iteratively until solving the auxiliary problem shows that there exists no variable with non-negative reduced cost. It is proven that the solution of the restricted master problem is now the optimal solution of the origin problem. [Rao 1969] proposes an integer programming formulation for cluster analysis by considering all possible clusters. With n objects, there are a total of 2 n possible clusters: C 1 to C 2 n . A matrix A = (a it ) with size n × 2 n is used to represent the relation between 4.3. Column generation framework for Constrained Sum of Squares Clustering 49

objects and clusters: a it = 1 if o i ∈ C t and a it = 0 otherwise. For the criterion of SSE, the cost c t for each cluster C t is dened as:

c t = SSE(C t ) = o i ∈Ct ||o i -m t || 2 = n i=1 ||o i -m t || 2 a it
where m t denotes the centroid of the cluster C t .

Let T be the set of all possible cluster indices, i.e. T = {1, 2, . . . , 2 n }. For all t ∈ T , let x t be boolean variables, x t = 1 denotes that the cluster C t is in the optimal partition.

The problem can now be formulated as an integer linear problem as:

min t∈T c t x t (4.1)
subject to:

             t∈T a it x t = 1 ∀i ∈ {1, . . . , n} (4.2a) t∈T x t = k (4.2b) x t ∈ {0, 1} (4.2c)
In this formulation, there are n constraints in Equation (4.2a) to express that, in the optimal partitions, each point must be in exactly one cluster. The constraint in Equation (4.2b) expresses that there are exactly k clusters. The dual of the master problem in (4.1) is expressed by [Rao 1969]:

max (-kσ + n i=1 λ i ) (4.3) subject to:            -σ + n i=1 a it λ i ≤ c t ∀t ∈ T (4.4a) λ i ≥ 0 ∀i = 1, 2, . . . , n (4.4b) σ ≥ 0 (4.4c)
The n dual values λ i correspond to n constraints in Equation (4.2a) and the dual value σ corresponds to the constraint in Equation (4.2b). Each column t in the master problem (4.1) corresponds to one constraint of the dual in Equation (4.4a).

This ILP formulation has 2 n variables, therefore it cannot be solved straightforwardly. The column generation method restricts the master problem in Equation The constraints can be among must-link, cannot-link and anti-monotone constraints. A constraint on a cluster is an anti-monotone constraint i whenever a cluster C satises this constraint, any sub-cluster C ⊆ C also satises it. According to the authors, an important observation to incorporate user-constraints to the column generation framework is: Mustlink, cannot-link and anti-monotone constraints can be evaluated on individual clusters in a partition. For example, a must-link constraint M L(o i , o j ) is satised in cluster C t if objects o i , o j are both in or not in this cluster. Similarly, a cannot-link constraint CL(o i , o j ) is satised in cluster C t if the cluster contains at most one of two objects. A must-link or cannot-link constraint is satised in a partition if it is satised in every cluster C t of the partition. With this observation, it is sucient to process user-constraints in the auxiliary solver: nding a column t with negative reduced cost such that C t satises all user-constraints. [START_REF] Babaki | [END_REF] proposes a BAB algorithm for nding the cluster C t which satises all user-constraints and corresponds to a negative reduced cost. The search process begins with an empty cluster C t and a set of candidate blocks such that covers all objects. A block is a set of objects that are connected by must-link constraints.

In case there is no must-link constraints, each block is equivalent to an object. The BAB algorithm performs a set-enumeration. In each iteration, a candidate block is considered to be aected to the cluster C t . With the notion of blocks, all must-link constraints are guaranteed to be satised in the returned cluster C t at the end. The cannot-link and anti-monotone constraints are used to remove incompatible candidate blocks.

Summary

In this chapter, we described recent declarative frameworks for solving data mining problems. We discussed the ideas and techniques involved and investigate advantages and disadvantages of proposed solutions. These researches inspire us that we could construct a general and exible framework for the problem of constrained clustering. Des améliorations du modèle sont étudiées, reposant sur des résultats théoriques, ce qui permet d'alléger le modèle sans changer la sémantique. De plus, parmi les critères d'optimisation, celui de la somme des dissimilarités se présente sous forme d'une somme linéaire. An d'avoir une propagation plus ecace pour ce critère, nous avons proposé un algorithme de ltrage qui exploite les informations données par une aectation partielle de variables.

Des expérimentations sur des bases de données classiques montrent la performance de notre modèle. De plus, la possibilité de combiner diérents types de contraintes utilisateur et la exibilité du choix du nombre de clusters permettent d'obtenir des solutions intéressantes.

Contents

Introduction

We introduce in this chapter our rst model for constrained clustering that can be provided to Constraint Programming solvers. We rst present the problem on which we focus. The problem consists of nding a partition composed of a given number of clusters that satises a set of user-constraints and optimizes a criterion. We dene the problem as follows:

Input

• O = {o 1 , . . . , o n }: set of points to be clustered. Points are indexed and named by their indices.

• k: number of clusters in a partition.

• A matrix of dissimilarities between two points. We denote by d ij the dissimilarity between points i and j.

Output A partition of O in k clusters C 1 , . . . , C k .
User-Constraints

∀o i ∈ C c , |{o j : o j ∈ C c , j = i, d ij ≤ } ∩ C c | ≥ M inP ts.
• Must-link constraint: A must-link constraint on two points i and j requires that these two points must be in the same cluster.

• Cannot-link constraint: A cannot-link constraint on two points i and j requires that these two points must be in dierent clusters.

Criterion

• No criterion: Finding all partitions with n points divided into k clusters that satisfy all given user-constraints.

• Diameter criterion: Finding a partition with n points divided into k clusters that satises all given user-constraints, and the partition has the minimal maximum diameter. The maximum diameter D of a partition is measured as:

D = max i,j s.t. λ i =λ j d ij
where λ i denotes the cluster that contains point i.

• Split criterion: Finding a partition with n points divided into k clusters that satises all given user-constraints, and the partition has the optimal minimum split. The minimum split S of a partition is measured as:

S = min i,j s.t. λ i =λ j d ij
• Within Cluster Sum of Dissimilarities (WCSD) criterion: Finding a partition with n points divided into k clusters that satises all given user-constraints, and the partition has the optimal WCSD. The WCSD is measured as:

W CSD = i,j s.t. λ i =λ j d ij
Let us notice that, for this criterion the dissimilarities are usually measured by the squared Euclidean distance.

For solving the described problem, we propose a model in Constraint Programming (CP). CP is a powerful paradigm for solving combinatorial optimization problems. The problems are encoded in CP by dening variables and constraints. A problem can be encoded in dierent ways in CP. CP relies on constraint propagation and search for solving the problem. As a result, the choice of variables and constraints is extremely important since dierent constraints have dierent ltering algorithms with dierent complexity and dierent levels of consistency. Finding a correct and ecient encoding in CP for a specic problem can be a dicult task.

CP provides a rich set of constraints that are designed to express relations between the variables. We present the constraints used in our model:

• Unary constraints: These constraints are dened on a single variable. The propagation of such a constraint is simple, as presented in Chapter 2.

• Binary constraints: Those are constraints between two variables x, y. In most CP solvers, binary constraints are guaranteed to be arc consistent with built in ltering algorithms. A basic arc consistency algorithm is presented in section 2.2 of Chapter 2.

• Global Cardinality Constraint: It bounds the number of times a value or a set of values is assigned to the variables in a set X to be at least l and at most u. We consider only a simple form of this constraint, which bounds the number of times a value v is assigned to variables in X. It is written as:

gcc(X, v, l, u)
where X is a set of variables, v is the value, l and u are the lower bound and upper bound of times v is assigned to variables in X. Another version used in this thesis is: gcc(X, y, l, u)

where y is a variable. This enforces that the value assigned to y must appear at least l and at most u times in the set X.

Example 17 Consider a set of variables X = {x 1 , x 2 , x 3 } with the domains Dom(x 1 ) = Dom(x 2 ) = Dom(x 3 ) = {1, . . . , 10}. The constraint gcc(X, 3, 1, 2) requires that at least one and at most two of the variables x 1 , x 2 , x 3 are assigned to the value 3.

Generalized arc consistency for this constraint can be achieved by a ltering algorithm, which is based on a network ow [Régin 1996]. We consider in the thesis two simple versions of the cardinality constraint: AtLeast and AtMost constraints, which are dened as:

AtLeast(X, v, m) AtM ost(X, v, m)
These constraints require that there are at least (at most) m variables in X that are assigned to value v.

• Element constraint: It gives access to a specied variable in an array of variables [x 1 , . . . , x n ]. It is written as:

element([x 1 , . . . , x n ], y, i)

It constrains the variable y to be equivalent to the i-th variable in the array, or: y =

x i . Generalized arc consistency for this constraint can be achieved [Gent et al. 2006].

• Reied constraint: It associates a boolean variable b to a constraint c. It is written as: b ↔ c (5.1)

In this constraint, the variable b takes the value 1 if the constraint c is satised and 0 otherwise. The propagation of reied constraints depends on the constraint c, and also depends on solvers. Most constraint solvers support reied versions of unary and binary constraints.

Example 18 Consider a set of variables X = {x 1 , x 2 , x 3 } with the domains Dom(x 1 ) = Dom(x 2 ) = Dom(x 3 ) = {1, . . . , 10}. Suppose that we need to post a constraint:

(x 1 = 2) → (x 2 = x 3 ) (5.2)
This constraint cannot be posted directly in most solvers. For modeling this constraint, we can use two auxiliary boolean variables b 1 and b 2 : Dom(b 1 ) = Dom(b 2 ) = {0, 1}. Then we post the 2 reied constraints as follows:

b 1 ↔ (x 1 = 2) (5.3) and b 2 ↔ (x 2 = x 3 ) (5.4)
Those are reied versions of binary constraints and are supported in most of constraint solvers. Finally, we post a binary constraint for binding b 1 and b 2 as follows:

b 1 → b 2 (5.5)
For instance, during the search, suppose that the domains of variables have been reduced to:

Dom(x 1 ) = 2, Dom(x 2 ) = {1, 3}, Dom(x 3 ) = {1, . . . , 10}
According to the constraint (5.3), the variable b 1 is instantiated: b 1 = 1. The propagation of constraint (5.5) then instantiates the variable b 2 : b 2 = 1. The propagation constraint (5.4) actives the constraint x 2 = x 3 , which restricts the domain of variable x 3 : Dom(x 3 ) = {1, 3}

Variables and domains

For modeling the problem described in section 5.1, we propose an encoding on two levels:

Representation of clusters For each cluster, a point in the cluster is chosen to be the representative point. The idea comes from the k-medoids algorithm [Kaufman & Rousseeuw 1987] where each cluster is represented by a medoid, i.e. the point whose average dissimilarity to all the objects in the cluster is minimal. However, identifying the medoid of each cluster requires to know the partition. We propose to replace the condition on representative points: a representative point is not necessary in the center of a cluster, but it is the point that has smallest index in the cluster. With this condition, modeling the representation of clusters in CP is simple, as presented by constraints given in the next section. This condition also helps to avoid symmetric solutions, as given a cluster, there is only one point satisfying the condition.

Representation of the assignment of points to a cluster Each cluster is associated with a representative point. We therefore propose to replace the assignment of points to a cluster to the assignment of points to the representative point. Suppose that the representative points are identied, the dissimilarity between a point and the representative of each cluster is known. Encoding the assignment of points to a representative point is intuitive and simple.

Based on these ideas, for modeling the problem of constrained clustering, we dene two sets of variables:

• I = {I 1 , . . . , I k }: For each cluster c ∈ {1, . . . , k}, the point with the smallest index is considered as the representative point of the cluster. The variable I c gives the representative point of the cluster c. The domain of each variable I c is all possible indices, hence Dom(I c ) = {1, . . . , n}.

• G = {G 1 , . . . , G n }: Assigning a point to a cluster becomes assigning the point to the representative of the cluster. For each point i ∈ {1, . . . , n}, the variable G i is the representative point of the cluster that contains the point i. At the beginning, any point can be a representative point, hence the domain of variable G i is also all possible indices, Dom(G i ) = {1, . . . , n}.

A complete assignment of those variables corresponds to a partition. We also propose another encoding which relies on only one one level, as presented in the next chapter. This partition corresponds to the assignment of variables:

I 1 = 1, I 2 = 3
(since 1 is the smallest index among {1, 2, 5, 6} and 3 is the smallest index among {3, 4, 7})

G 1 = G 2 = G 5 = G 6 = 1
(since 1 is the representative of the rst cluster)

G 3 = G 4 = G 7 = 3
(since 3 is the representative of the second cluster). A oat variable is introduced for representing the optimization criterion. It is denoted by D for the maximal diameter, S for the minimal margin and V for the Within Cluster Sum of Dissimilarities. The domains of D and S are the interval whose lower (upper) bound is the minimal (maximal, resp.) dissimilarity between any two points. The domain of V is upper-bounded by the sum of the dissimilarities between all pairs of points.

Constraints

Modeling Partition Constraints

The relationships between points and their clusters are expressed by the following constraints:

The representative of a representative point is itself Obviously, if a point i is the representative point of the cluster c, then I c = i. The variable G i must be assigned to i (G i = i) because the point i is in the cluster c, and the point i is also the representative of this cluster. As a result, the variables G Ic , giving the cluster to which the representative point I c belongs to, must be equal to I c . This can be expressed by an element constraint. For each value c ∈ {1, . . . , k}, we post a constraint:

element([G 1 , . . . , G n ], I c , I c ) (5.6)
With those constraints, as soon as a representative point i of a cluster c is identied, both variable I c and G i are instantiated to the value i:

I c = G i = i.
The representative of a cluster must be the smallest index in the cluster This condition makes the representative be unique for each cluster. For each point i, the value G i giving the representative point of its cluster, must be smaller or equal than i. To represent this constraint, for each value i ∈ {1, . . . , n}, we post a unary constraint:

G i ≤ i. (5.7)
Example 20 Given the dataset in Example 19, the constraints (5.7) restrict the domain of variables in G to:

Dom(G 1 ) = {1}, Dom(G 2 ) = {1, 2}, . . . , Dom(G 7 ) = {1, . . . 7}.
The representative of each point must be one and only one of all the representatives The value of a variable G i must be among the possible representative points. It means that G i must be assigned to the value of exactly one of the variables in I. This relation can be expressed by a global cardinality constraint. For each value i ∈ {1, . . . , n}, we post a constraint: gcc(I, G i , 1, 1) (5.8) This means that the value of G i must appear exactly once time in I. With these constraints, when all variables in I are assigned, the domains of variables in G contain only the values of representative points.

Example 21 Given the dataset in Example 19, suppose that the two representatives points are identied: I 1 = 1, I 3 = 3. With the constraints (5.6), (5.7) and ( 5.8) the domains of all variables in G are ltered to:

Dom(G 1 ) = Dom(G 2 ) = {1}, Dom(G 3 ) = {3}, Dom(G 4 ) = Dom(G 5 ) = Dom(G 6 ) = Dom(G 7 ) = {1, 3}.
Symmetry breaking constraints In our model, given an assignment of variables in I and G, any permutation of the values of the representative points leads to the same partition. Figure 5.2 presents an example of symmetric solutions, if there is no restriction on the choice of representative points. A simple way to avoid symmetry solutions is to impose a strict ordering on the representatives: the representatives must be in an ascending order. This relation can be expressed by binary constraints. For each value c ∈ {1, . . . , k -1}, we post a constraint:

I c < I c+1 .
Since a representative must have the smallest index, the representative of the rst cluster must be the rst point, so we post a constraint:

I 1 = 1.

Modeling the user-constraints

The advantage of modeling with CP is that many dierent types of user constraints can be directly added. Both cluster-level constraints and instance-level constraints can be formulated within our model. All popular user-dened constraints may be straightforwardly integrated: This constraint expresses that there must be at least α variables in G taking the value of the variable I c , which is equivalent to there are at least α points in the cluster c. This constraint requires that there must be at most β variables in G taking the value of the variable I c .

Maximal size β of clusters

Split δ-constraint It expresses that the split between any two clusters must be at least δ. This constraint is satised if for any two points i and j, if the dissimilarity between them is less than δ, we require them to be in the same cluster. Therefore, for each i < j ∈ [1, n] satisfying the condition d ij < δ, we put a binary constraint:

G i = G j .
Diameter constraint γ It expresses that the diameter of each cluster must be at most γ. This constraint is satised if for any two points i, j, if the dissimilarity between them is greater than γ, we require them to be in dierent clusters. Therefore for each i < j ∈ [1, n] such that d ij > γ, we put a binary constraint:

G i = G j
Density constraint It requires that for each point i, its neighbourhood of radius contains at least M inP ts points belonging to the same cluster as i. To express this constraint in our model, for each point i ∈ [1, n], we compute the set of variables X i composed of variables G j (j ∈ {1, . . . , n}) such that d ij < . We then post an AtLeast constraint: AtLeast(X i , G i , M inP ts).

It requires that, in the set of variables X i , there must be at least M inP ts variables that have the same value as G i , which means that there must be at least M inP ts points in the same cluster as point i.

Must-link constraint A must-link constraint on two points i and j is expressed by a binary constraint:

G i = G j .
Cannot-link constraint A cannot-link constraint on i and j is expressed by a binary constraint:

G i = G j .

Criterion Modeling

As presented in Chapter 2, a Constraint Optimization Problem (COP) is a CSP together with an objective function f () to be optimized. Typically, a COP is solved in CP by adding a constraint whenever a solution ∆ is found during the search. The constraint requires that the value of the optimization function must be better than f (∆) in any future solution. For modeling a COP, in general, a variable, i.e. x, is dened and users have to post constraints on x to express the function f (). Then an optimization constraint is posted, e.g. in case of minimization:

minimize x
This constraint ensures that, whenever a solution ∆ is found, a new constraint is added to the model, e.g. in case of minimization:

x < f (∆) (5.9) where f (∆) is the actual value of the function f () on the solution. f (∆) is obtained by the value of the variable x in the solution. In some problems, the variable x is not instantiated and f (∆) is obtained by a lower or upper bound of x. By using the constraint (5.9), the last found solution is guaranteed to be the solution with the minimal value of the function f () among all found solutions. When a complete search is used, the last found solution is guaranteed to be the optimal solution. Note that the constraint (5.9) is strict, providing an increasingly tight bound. As a result, only one optimal solution is found, even if there exists dierent optimal solutions with the same value of f ().

Modeling the diameter criterion

The criterion of diameter is encoded with the oat variable D. At the beginning, the domain of D is not a set of values, but an interval:

Dom(D) = [d min , d max ]
5.4. Criterion Modeling 61 where d min and d max are the minimum and the maximum dissimilarity between any two points. To express the criterion, a straightforward way to encode it is to put a requirement on each pair of points: if two points i and j are in the same cluster, the dissimilarity between these two points should be equal or inferior to the value of the maximum diameter:

G i = G j → d ij ≤ D (5.10)
With these constraints, whenever two points i and j are known to be in the same cluster, the lower bound of the variable D is modied: D.lb = max(D.lb, d ij ). When all variables in G are instantiated, the lower bound of D is the maximum diameter.

Note that the constraint (5.10) is a logic constraint in the form a → b, which is equivalent with !b →!a, with the meaning: if the dissimilarity between any two points is superior to the diameter, then they must be in dierent clusters.

d ij > D → (G i = G j ) (5.11)
The condition d ij > D is satised if the dissimilarity d ij is superior to the upper bound of the variable D. Whenever d ij > D.ub, the constraint G i = G j has to be activated to restrict the domain of variables G i and G j .

For modeling the constraints (5.10, 5.11), we propose to use reied versions of binary constraints. For each pair of points i and j, two auxiliary boolean variable b ij1 and b ij2 are generated and three constraints are posted. The rst one models the condition that the dissimilarity between any two points is superior to the diameter:

b ij1 ↔ d ij > D.
(5.12)

The second one models the condition that the two points are not in the same cluster:

b ij2 ↔ (G i = G j ).
(5.13)

The third one binds the two conditions:

b ij1 → b ij2 . (5.14)
The propagation of the constraints (5.12, 5.13, 5.14) is as follows:

• If b ij1 = 1, which happens when d ij > D.ub holds in the constraint (5.12), then the variable b ij2 is instantiated to 1 because of the constraint (5.14). The reied constraint (5.13) will then propagate the binary constraint

G i = G j . • If b ij1 = 0 or b ij1 is not instantiated (b ij1 ∈ [0, 1]
), there is no propagation.

• If b ij2 = 1 or (G i = G j ) holds, there is no propagation.

• If b ij2 = 0 or (G i = G j ) holds in the constraint (5.13), then the variable b ij1 is instantiated to 0 because of the constraint (5.14). With the reied constraint (5.12), the constraint d ij ≤ D is propagated. This is done by setting the lower bound of the domain of D: D.lb = max(D.lb, d ij )

These propagations ensure that the lower bound of D is always the maximum diameter of clusters. Since G 2 = G 3 , we can see that b 24 and b 34 must be equal and then must not be equal to 1, otherwise the constraint is violated. We should then infer that point 4 cannot be in the cluster number 2, as points 2 and 3, that means value 2 should be removed from the domain of G 4 . This ltering however is not done, since the constraints are considered independently.

Several recent works have proposed more ecient ltering for the sum constraint, when it is considered with other constraints. For a sum constraint y =

x i with inequality constraints x j -x i ≤ c, a domain ltering algorithm reduces the domain of x i when new bounds for y are known [START_REF] Régin | [END_REF]. A bound-consistency algorithm is proposed for a sum constraint with increasing order constraints x i ≤ x i+1 [Petit et al. 2011] or with a constraint alldierent(x 1 , . . . , x n ) [Beldiceanu et al. 2012]. These cases however do not t the WCSD criterion constraint (5.22). A generic bound-consistency algorithm for a sum constraint with a set of constraints is proposed in [Régin & Petit 2011]. In our case, the domain of G i is a set of representative indices, which is not an interval in general, and from which we wish to remove any inconsistent value. We have developed a ltering algorithm for a new global constraint on a variable V , an array of variables G of size n and an array of constants d, which is of the form:

V = i,j∈[1,n],i<j (G i = G j )d ij .
(5.23)

Taking into account the partitioning problem, the domain of each variable G i is the set of the representative indices of the clusters that can contain point i. Let us assume that the domain of the variable V is [V.lb, V.ub) where V.lb is the lower bound, which can be initially 0, and V.ub is the upper bound, which can be the value of V in the last solution with a branch-and-bound search. Suppose that we have a partial assignment of variables in G, where there is at least one point assigned for each group (e.g the representative of the group). Let K be the set of points i which have been already assigned to a group (G i is instantiated) and U the set of the unassigned points. The sum in (5.23) is split into three parts

V = V 1 + V 2 + V 3
, where:

• V 1 is the sum of dissimilarities between the assigned points:

V 1 = i,j∈K,i<j,G i =G j d ij
• V 2 is the sum of dissimilarities between the unassigned points and the assigned points:

V 2 = i∈U,j∈K (G i = G j )d ij
• V 3 is the sum of dissimilarities between the unassigned points:

V 3 = i<j,i,j∈U (G i = G j )d ij
The evaluation of these three parts is based on [Klein & Aronson 1991]. The value of V 1 can be calculated exactly because the set K is already known. For the second part, the value of V 2 is unknown because of the unassigned points. However, a lower bound of V 2 , denoted by V 2 .lb, can be calculated by a sum of the minimum contribution of all the unassigned points. Since each unassigned point i will be assigned to a group, it will contribute to V by a sum of dissimilarities between point i and all the points of K that are in that group. The minimal contribution v 2i of the point i is the minimal amount when considering all k groups, with respect to the assigned points:

v 2i = min c∈[1,k] ( j∈K∩Cc d ij ).
A lower bound of V 2 is then the sum of v 2i :

V 2 .lb = i∈U v 2i .
For the third part, the value of V 3 is unknown too, but a lower bound of V 3 can be calculated by a heuristic. We recall that V 3 is the sum of all d ij with i and j in the same group. Let p = |U |, the minimal number of terms d ij in the sum V 3 is the minimal number of within-group pairwise connections 1 , while considering all partitions of p points into k groups. This number can be calculated by a function, i.e. f (p, k). With a set U of unassigned points, with the constants d ij (i, j ∈ U ) ordered increasingly, a lower bound V 3 , denoted by V 3 .lb, is then calculated by the sum of the f (|U |, k) rst constants in this order.

Example 22 With p = 10, k = 3 and with a partition into 3 groups of sizes 2, 3 and 5, the number of within-group pairwise connections is 14. The minimal value of this number is 12, corresponding to a partition into 3 groups of sizes 3, 3 and 4.

Theorem 22 Let m be the quotient of the division of p by k and m the remainder. The minimal total number of connections for all groups can be calculated as follows:

f (p, k) = (km 2 + 2mm -km)/2.
Proof Let the number of points in each group c be m + α c , with α c < 0 when the group c has less than m points, α c ≥ 0 otherwise. We have then The number of pairwise connections in a group c is:

1≤c≤k (m + α c ) = p = km + m .
(m + α c )(m + α c -1)/2.
The total number for all groups is:

1≤c≤k (m + α c )(m + α c -1)/2 = ( 1≤c≤k (m + α c ) 2 -1≤c≤k (m + α c ))/2 = (km 2 + 2mm + 1≤c≤k (α 2 c -km -m ))/2
By replacing the value of m (m = 1≤c≤k α c .), we have:

m ≤ 1≤c≤k |α c | ≤ 1≤c≤k α 2 c
where α c are integers. Therefore the total number for all groups is greater or equal to:

(km 2 + 2mm -km)/2.
The equality is reached when α c is 1 for m groups and is 0 for k -m groups

Example 23 Given a dataset with 8 points as illustrated in Figure 5.4 and a clustering task dividing the points into 2 groups (green and red). Suppose that there are 5 assigned points (3 red points and 2 green points) and 3 unassigned points. The value of V 1 is calculated exactly by the sum of the solid black lines. The lower bound V 2 .lb is the sum of the dash red lines and dash green lines. With 3 unassigned points, we have p = 3, k = 2, m = 1 and m = 1, the minimum total number of connections is

f (p, k) = f (3, 2) = (km 2 + 2mm -km)/2 = 1.
Therefore, the lower bound V 3 .lb is the dot blue line.

A lower bound of variable V is given by:

V.lb = V 1 + V 2 .lb + V 3 .lb
This lower bound is used for two purposes:

• Detecting the failure during the branch-and-bound search, it happens when V.lb ≥ V.ub.

• Filtering inconsistent values of unassigned variables. For each value v of an unassigned variable G i , a new lower bound, denoted by V .lb, is computed with the assumption G i = v (the computation is explained in the next paragraph). This value is inconsistent if V .lb ≥ V.ub. For example in Figure 5.3, value 2 can be removed from the domain of variable G 4 because with the assumption G 4 = 2, the new lower bound V .lb = 6, and it is greater than the upper bound V.ub = 5.

V.lb = V 1 + V 2 .lb + V 3 .lb
The ltering algorithm is presented in Algorithm 11. This algorithm uses two arrays add and min:

• add[i, c] is the added amount if i is assigned to group c:

add[i, c] = j∈K∩Cc d ij .
• m[i] is the minimal added amount while considering all possible assignments for i:

m[i] = min c add[i, c].
Since the constants d ij must be ordered increasingly in the computation of V 3 .lb, they are ordered once in the array ord, so ord[pos] gives the constant d ij in the order at position pos, and px[pos] (py[pos]) gives the index i (j, resp.) of the constant. For the time being, the ltering algorithm is developed for the clustering task, where values in the domain of G i are the representatives of all clusters for which point i can be assigned. Given a value v in the domain of a variable G i , the function gr(v) gives the index of the cluster corresponding to v.

The lower bound of V is revised in line 26. Line 30 lters the domain of G i (i ∈ U ): for each value v in the domain, in case of assignment of i into group gr(v), a new lower bound for V.lb is V .lb = V 1 + V 2 .lb + V 3 .lb with: add[i, gr(v)] because point i is supposed to be assigned to group gr(v), the sum of dissimilarities between instantiated points increases of add [i, gr(v)].

• V 1 = V 1 +
• V 2 = V 2 .lb -m[i]
because point i is no more unassigned, the contribution of point i in the calculation of V 2 .lb must be removed.

• V 3 .lb is the sum of the rst f (|U | -1, k) elements of ord that are related to U \ {i}.

In order to reduce the complexity of the ltering algorithm, we use V 4 instead of V 3 .lb. Here, V 4 is the sum of the rst f (|U | -1, k) elements of ord (possibly some related to i). It is obvious that

V 3 .lb ≥ V 4 .
The new lower bound is: is considered for branching. This point is in the frontier between two clusters. Although this point is closer to the cluster on the left, there is a high probability that it is not true in the optimal solution, as the dierence between the dissimilarities of this point to the cluster on the left and on the right is close to zero. If assigning the point B to the cluster on the left is not correct, we want to detect this failure soon. Branching on point B makes the rst part (V 1 ) of the calculation on bounds on WCSD increase a lot, and the lower bound of W CSD is larger than the bound when we make the branching on point A. As a result, the ltering algorithm for WCSD is more eective, we have a better chance to detect a failure and to lter more values in the domain of variables in G.

(V 1 + add[i, gr(v)]) + (V 2 .lb -m[i]) + V 3 .lb

Model Improvements

By using a complete branch and bound search, our model allows to nd the optimal solution. In order to improve the eciency, dierent aspects have been studied.

Search strategy improvement by point ordering

In order to instantiate cluster representatives, the variables in I are taken in the order I 1 , . . . , I k when enumerating with values. Since a representative must have the smallest index in the cluster, for each variable, values (point indexes) are enumerated in the as- We present in Algorithm 12 the reordering algorithm. It is based on the FPF heuristic [Gonzalez 1985] to reorder points. The idea is to use the FPF algorithm with k = n (each point is in a cluster), so that each point will be taken in turn. The order in which the algorithm chooses points gives the order of points. Figure 5.7 presents points indices before and after the reordering. Figure 5.8 presents the rst solution found with the old and new ordering. With the new ordering, the rst solution usually has a diameter which is close to the optimal value. It is very important as this value is used for the upper bound of the variable D.

Constraint improvements

For a given k, without any user-constraints, the diameter d F P F returned by the FPF algorithm has been proved by [Gonzalez 1985] to satisfy:

d opt ≤ d F P F ≤ 2d opt (5.24)
where d opt is the optimal diameter, which is aimed to be found. This entails bounds for the variable D, whose domain is

[d F P F /2, d F P F ].
Moreover, if the dissimilarity between any two points i, j is greater than d F P F , we can directly add constraints to indicate that they must be in dierent clusters. It replaces the constraints

d ij > D → (G i = G j ).
That is, for all points i < j ∈ [1, n]:

• If d ij < d F P F /2:
There is no need to post the constraints (5.12, 5.13, 5.14) for the criterion of diameter. The reason is that as d ij < d F P F /2 and d F P F /2 ≤ d opt , the condition d ij > D is always false, there is no propagation on these constraints even if we post them.

Algorithm 12: Reordering algorithm based on the FPF heuristic Input: a number of points n, a number of clusters k, a matrix of dissimilarity

(d ij ) Output: a new index of points in array N ewIndex[1..n] i ← the furthest point (i ← arg max j t∈[1,n] d jt ) point[i] ← visited N ewIndex[i] ← 1 foreach j ∈ [1, n] do if point[j] is not visited then DistanceT oV isitedP oint[j] ← d ij end end foreach index ∈ [2, n] do
i ← the point that is not visited and has the largest value DistanceT oV isitedP oint

[i] point[i] ← visited N ewIndex[i] ← index foreach j ∈ [1, n] do if point[j] is not visited then DistanceT oV isitedP oint[j] ← d ij end end end • If d F P F /2 ≤ d ij ≤ d F P F :
There is no change, three constraints (5.12, 5.13, 5.14) are posted for the criterion of diameter.

• If d ij > d F P F : a constraint G j = G i is posted, instead of three constraints (5.12, 5.13, 5.14), because the condition d ij > D is always true.

Lots of reied constraints can be removed from the model, which improves the model eciency.

In case there are user-constraints, the properties: d opt ≤ d F P F is no longer true. The reason is that the optimal diameter in case there are user-constraints is actually equal or superior to the optimal diameter in case without user-constraints. However, we always have:

d F P F ≤ 2d opt or d F P F /2 ≤ d opt .
As a result, only the lower bound of D is changed to d F P F /2. For all points i < j ∈ [1, n], we make the changes:

• If d ij < d F P F /2:
There is no need to post the constraints (5.12, 5.13, 5.14) for the criterion of diameter.

• If d F P F /2 ≤ d ij : Three constraints (5.12, 5.13, 5.14) are posted for the criterion of diameter.

Experiments

Our model is implemented on Gecode solver version 4.2.1. Twelve databases from the repository UCI [START_REF] Bache | [END_REF] are used in our experiments. They vary by their size, number of attributes and number of classes. Table 5.1 summarizes information on these datasets, presented in the increasing order of the number of objects. The experiments are all performed on a 3.4GHz Intel Core i5 processor running Ubuntu 12.04. We experiment the performance of our model with each criterion by comparing with other exact methods. To our knowledge, there is no exact algorithm for these criteria that supports any kind of user constraints for a general value k ≥ 3. For the criterion of diameter, our model is compared with the branch-and-bound approach [Brusco & Stahl 2005] and the algorithm based on graph coloring [Delattre & Hansen 1980], in the case without userconstraints. For the split criterion, nding a partition that optimizes the criterion is a polynomial problem. We experiment our model adding user-constraints on the diameter.

For the criterion of WCSD, our model is compared with the Repetitive Branch-and-Bound Algorithm (RBBA) [Brusco & Stahl 2005]. These algorithms are detailed in section 3.3.2 of Chapter 3. We also demonstrate the benet of using user-constraints in our model.

Minimizing the maximal diameter of clusters

We compare our previous model (denoted by CP) with the branch-and-bound approach [Brusco & Stahl 2005] (denoted by BaB) and the algorithm based on graph coloring [Delattre & Hansen 1980] (denoted by GC). Our programs are available at httpXGGpRlusteringFom. The program BaB can be found at httpXGGmilerFfsuFeduG£mrusoG whereas the program GC is not available and we have coded it in C++ using a well known available graph coloring program [START_REF] Mehrotra | [END_REF]. We consider clustering without user-constraints since the other 5.2: Performance with the criterion of minimizing the maximal diameter algorithms cannot handle them and in the best of our knowledge, there is no exact algorithm handling user-constraints with this criterion. In the experiments, the time-out is set to 1 hour and the Euclidean distance is used to compute the dissimilarity between objects. The number of clusters k is set to the number of real classes given in Table 5.1. Table 5.2 shows the results of our experiments. For each dataset we present the value of D opt (the optimal diameter) in the second column and the run time in seconds of each system. Thesign is used when the system cannot complete the search after 1 hour. All the algorithms are exact and therefore all nd the same value for the optimum diameter.

It is clear that with these datasets, our model is the most ecient in all cases and can nd the optimal diameter in 10 datasets. Among the programs, BaB algorithm is the least ecient, it is not able to solve datasets with more than 300 objects. The performance of GC is better than that of BaB but it decreases rapidly if the number of objects n is over 500. The BaB algorithm is based on the bounds of the maximal diameter to detect

- - Multi Features - - Image Segmentation - - Waveform - -
Table 5.3: Maximizing the minimal split between clusters with a diameter constraint failures during the search, while the GC algorithm considers all available dissimilarities in descending order to nd the optimum diameter. Our model, which exploits the benets of Constraint Programming such as constraint propagation and appropriate search strategies, is more ecient.

Maximizing the minimal split between clusters

Finding a partition maximizing the split between clusters is a polynomial problem. However, this is no longer true with user constraints, as for instance with a diameter constraint.

To our knowledge, there is no exact algorithm for this criterion that supports any kind of user constraints for a general value k ≥ 3. We experimented our model by adding a diameter constraint. In order to set it we use the results given in Table 5.2: the diameter of each cluster must not exceed 1.5D opt . The number of clusters k is set to the actual number of classes given in Table 5.1. The results are given in Table 5.3. For each dataset, we present the optimal split S opt and the total execution time in seconds. Our model is able to solve datasets with less than 1000 objects.

Minimizing the Within Cluster Sum of Dissimilarities (WCSD)

Finding an exact solution for minimizing the WCSD is dicult and we compare our model with the Repetitive Branch-and-Bound Algorithm (RBBA) [Brusco & Stahl 2005]. The program is available at httpXGGmilerFfsuFeduG£mrusoG, which, to our best knowledge, is the best exact algorithm for the WCSD criterion. In these experiments, the dissimilarity is measured as the squared Euclidean distance.

Minimizing the WCSD is a dicult task since the propagation of the sum constraint is less ecient than the propagation of the diameter constraint. Without user-constraints, both our model and the RBBA approach can nd the optimal solutions only with the Iris dataset. Our model needs 4174s to complete the search whereas RBBA takes 3249s. However, with appropriate user-constraints, the performance of our model can be signicantly improved. WCSD and the split constraint Let us add a split constraint δ, where δ ranges from 0% (no constraint) to 12% of d max , where d max is the maximum dissimilarity between two objects in the dataset. Table 5.4 (left) reports the WCSD value of an optimal solution and the total computation time. It shows that when the split constraint is weak, the optimal WCSD value does not change but the computation time decreases signicantly when this additional constraint becomes stronger. The reason is that the total number of feasible solutions decreases and the search space is reduced. When the split constraint is weak, propagating this constraint is more time-consuming than its benets.

WCSD and must-link constraints Let us now add must-link constraints, where the number of must-link constraints, generated from the true classes of objects, varies from 0.2 to 1% of the total number of pairs. The results are expressed in Table 5.4 (right), giving the WCSD value and the total computation time. In fact, the optimal value of WCSD, with no information on classes, does not correspond to the WCSD found when considering the partition of this dataset into the 3 dened classes. The more must-link constraints, the less computation time is needed for nding the optimal value, and the closer to the value of WCSD, when considering the 3 initial classes. The reduction of computation time can be easily explained, since when an object is instantiated, objects that must be linked to it are immediately instantiated too. Furthermore, with any kind of additional constraint, the total number of feasible solutions is always equal or less than the case without constraint.

WCSD and appropriate user-constraints Finding an exact solution minimizing the WCSD is dicult. However, with appropriate combination of user-constraints, the performance can be boosted. Table 5.5 presents some examples where our model can get an optimal solution with dierent user-constraints, which reduce signicantly the search space.

Flexibility of model

Our system nds an optimal solution when there exists one; otherwise our system proves that there is no solution. Let us show the interest of combining dierent kinds of constraints. Figure 5.9 presents 3 datasets in 2 dimensions, similar to those used in [Ester et al. 1996]. The rst dataset is composed of 4 groups with dierent diameters. The second one is more dicult, since groups do not have the same shape. The third one contains outliers: outliers are not handled and are therefore integrated in classes.

When optimizing the maximal diameter, the solver tends to nd rather homogeneous groups, as shown in Figure 5.10. Adding a split constraint (with δ = 5% of the maximum dissimilarity between pairs of points) improves the quality of the solution (see Figure 5.11). Let us notice that maximizing the minimum split allows also to nd this solution.

Concerning the third dataset, minimizing the maximum diameter or maximizing the minimum margin do not allow nding a good solution (see Figure 5.12). The quality of the solution is improved when a density constraint is added with M intP ts = 4 and = 25% from the maximal dissimilarity between pairs of points.

Summary

Modeling a problem in Constraint Programming in an ecient and eective way requires a lot of skills and eorts. In practice, in order to nd a good model, we may have to try many models before achieving a good formulation. In this chapter we have proposed a general and declarative model in Constraint Programming for solving the problem of Constrained Clustering. Our model allows to choose among dierent optimization criteria and to integrate various kinds of user-constraints. Experiments on various classical datasets show that our model is competitive with other exact approaches. To increase the eciency and to deal with larger datasets, we address some analysis that may help us to improve the model:

• The model requires to specify the number of clusters k, which is not always known in the task of clustering. In order to allow more exible models, where the number of clusters is not xed, we may need other representations of variables.

• For the criterion of maximum diameter and minimum split, we have to post a large number of reied constraints. Many of these constraints may never help the solver to lter inconsistent values in the domain of variables. For example, for the criterion of minimum split, if the dissimilarity between two points i and j is large, it is obvious that this dissimilarity has nearly zero chance to be the optimal split. As the result, the reied constraints on these points cannot remove values in the domain of variables. However the constraint solver still has to check these constraints every time it visits a new node in the search tree. Finding other ways for modeling the criteria may speed up the performance.

Introduction

In this chapter, we propose a second model, still based on Constraint Programming, but signicantly dierent from the one presented in Chapter 5. In the previous model, in order to identify a partition, two sets of variables were introduced, namely a variable for each cluster identifying a cluster by one of its points and a variable for each point assigning it to a cluster. The number of classes had to be xed beforehand. The new model we present here contains only a variable for each point, giving the number of the cluster the point belongs to. The model is thus lighter in terms of the number of variables, it also enables to remove the restriction on the number of classes, only bounds on the number of classes are required. We develop dedicated ltering algorithms for global constraints for the criteria of minimizing the maximal diameter and maximizing the minimal split between clusters.

Variables and domains

We present our second CP model, which is based only on a set of variables for the assignment of a number of class to each point. In this model, the number of classes k is only bounded by k min and k max , where k min and k max are given by the user. The dierence is signicant, the new model oers more exibility, while it is more ecient than the previous model.

In this model, the clusters are identied by their index, which varies from 1 to k for a partition into k clusters. To represent the assignment of points to clusters, we use a set of integer variables:

G = {G 1 , . . . , G n }
The domain of each variable G i is all possible cluster indices, which ranges from 1 to k max :

Dom(G 1 ) = • • • = Dom(G n ) = {1, . . . , k max }
An assignment G i = c means the point i is put into the class number c. A complete assignment of those decision variables corresponds to a partition. For modeling the criteria, we use the same variables as in the previous model: D for the maximal diameter criterion, S for the minimal split criterion and V for the Within-Cluster Sum of Dissimilarities criterion. Figure 6.1 presents a sample data and an assignment of variables in G which corresponds to a partition. The left image presents the initial domains of the variables in G in which we want to solve a clustering task into 1 or 2 clusters. The right image presents an assignment that corresponds to a partition of 2 clusters.

Constraints

Modeling Partition Constraints

Similarly to the previous model, we need constraints to ensure that each possible partition with at least k min clusters and at most k max cluster corresponds to only one complete assignment of variables in G and vice versa.

Without any constraint, given a partition with k clusters (k min ≤ k ≤ k max ), there exists dierent assignments of variables which correspond to the same partition.

Dom(G 1 ) = . . . = Dom(G 7 ) = {1, 2, 3}
Let us consider a partition of 2 clusters, the rst one composed of o 1 , o 2 , o 5 and o 6 and the second one composed of the remaining points. There are in total 6 assignments of variables, presented in Figure 6.2, which express the same partitions.

From an assignment, we can obtain a symmetric solution by:

• Permutation of the values of variables, for example the solution A and B, C and D, E and F in Figure 6.2.

• Changing the values of variables corresponding to a cluster to another cluster index value that has not yet been assigned, for example the solution A and C, C and E in Figure 6.2.

To avoid the symmetric solutions C, D, E and F , it is simple to put a constraint that: if there are k clusters, then for each value c ∈ {1, . . . , k}, there exists at least one variable G i that is assigned to c. To avoid the symmetry solution B, a simple way is to post a condition that: For any two cluster indices c < c that appear in an assignment, if a variable G i is assigned to c , there must exists a variable G j that is assigned to c and j < i.

The global constraint precede [Law & Lee 2004] helps to achieve both requirements above. For avoiding all possible symmetry solutions, all we need is to post a constraint:

precede([G 1 , . . . , G n ], [1..k max ]).
This constraint imposes that G 1 = 1 and moreover, if G i = c with 1 < c ≤ k max , there must exist at least an index j < i with G j = c -1.

The fact that there are at least k min clusters means that all the numbers among 1 and k min must be used in the assignment of the variables G i . When using the constraint precede, one only needs to require that at least one variable G i is equal to k min . This is expressed by an AtLeast constraint:

AtLeast(G, k min , 1)
Since the domain of each variable G i is [1, k max ], there will be at most k max clusters. If the user needs exactly k clusters, all he/she has to do is to set k min = k max = k.

Modeling the user-constraints

In our previous model, the user-constraints can be straightforwardly integrated with constraints on variables in G and I. Although in this new model, the variables in G express the cluster indices, not the representative, the underlying meaning is the same in the two models: two points i and j are in the same cluster if and only if G i = G j . There are therefore not many dierences in the modeling of user-constraints.

Minimal size α of clusters It requires that each cluster has a number of points greater than a given threshold α. This constraint can be expressed by a constraint AtLeast. If the number of clusters is xed to k, the constraint can be simply expressed by global cardinality constraints: for each value c ∈ {1, . . . , k}, we post a constraint:

AtLeast(G, c, α).
However, in our model, the number of clusters is bounded by two values k min and k max . We, therefore, propose to express the minimal size constraint by n constraints AtLeast: for each value i ∈ {1, . . . , n}, we post a constraint:

AtLeast(G, G i , α).
Each constraint expresses that there must be at least α variables in G taking the value that is assigned to G i . This is equivalent with there are at least α points in the cluster that contains point i. This constraint requires that there must be at most β variables in G taking the value c.

Maximal size β of clusters

Split δ-constraint It expresses that the split between any two clusters must be at least δ. Similar to the previous model, for each i < j ∈ [1, n] satisfying the condition d ij < δ, we put a binary constraint:

G i = G j .
Diameter constraint γ It expresses that the diameter of each cluster must be at most γ. For each i < j ∈ [1, n] such that d ij > γ, we put a binary constraint:

G i = G j .
Density constraint It requires that for each point i, its neighbourhood of radius contains at least M inP ts points belonging to the same cluster as i. Similar to the previous model, for each point i ∈ [1, n], we compute the set of variables X i which contains variables G j of points j ∈ {1, . . . , n} such that d ij < . We then post a constraint AtLeast:

AtLeast(X i , G i , M inP ts).

It requires that, in the set of variables X i , there must be at least M inP ts variables that is equal to G i , which means there must be at least M inP ts points in the same cluster with point i.

Must-link constraint A must-link constraint on two points i and j is expressed by a binary constraint:

G i = G j .
Cannot-link constraint A cannot-link constraint on i and j is expressed by a binary constraint:

G i = G j .
6.4 Criterion Modeling

Modeling the diameter criterion

In our previous model, the criterion of diameter is expressed by implication constraints:

d ij > D → (G i = G j ) (6.1)
For modeling each implication constraint, in most constraint solvers, we have to use several primitive constraints. For instance, we proposed to use two reied constraints and one binary constraint for modeling each implication constraint, as presented in the previous chapter. In order to have better interactions and propagation between these relations, a better way is to sum up these relations into one global constraint with a dedicated ltering algorithm.

We have developed a global constraint for the criterion of maximum diameter, de-

Modeling the split criterion

For the criterion of split, we have developed a similar global constraint split(G, S, d) ensuring that for every couple 1 ≤ i < j ≤ n:

S > d ij → G i = G j . (6.2)
The ltering algorithm is presented in Algorithm 14, where Dom(S) = (S.min, S.max].

In this algorithm, if S.min has been changed and if S.min ≥ d ij , the relation (6.2) will infer G i = G j , which is propagated by enforcing Dom(G i ) = Dom(G j ). Otherwise, in line 18, in case some variables G j have been instantiated, if G i = G j by ( 6.2) we infer S ≤ d ij , so the upper bound of S can be changed. The worst case complexity is O(n 2 ).

assignment of a point to a new cluster is favoured if there are unused cluster numbers. In this case, the least remaining number is chosen. The closest cluster c to the point i is chosen and two alternatives are created G i = c and G i = c. This strategy, which takes into account all the points already assigned to decide the branching, is dierent from the one used in the previous model, where the branching depends on the dissimilarity between the point i and the representative of each cluster.

For the criterion WCSD The same strategy is used as in the previous model. For nding the rst partition, the chosen variable G i and value c are those such that the assignment G i = c increases WCSD as least as possible. After nding the rst solution, the strategy is changed following the rst fail principle. A value s ic for each point i and each cluster c is dened as the sum of the squared dissimilarity between i and all points j already assigned to the cluster c. At each branching point, for all points i with G i unassigned, the minimal value s i = min c∈Dom(G i ) s ic is computed. The variable G i with the smallest value s i is then chosen and the value c = arg min s ic is chosen. Two alternatives where G i = c and G i = c are created for the branching.

6.6 Experiments

Comparison of two models

We compare our rst model (denoted by CP1) with our second model (denoted by CP2) in order to study the following questions:

• Q1: Do both models eliminate all symmetry solutions?

• Q2: Which one better models variables and partition constraints?

• Q3: How do the dedicated global constraints for the criteria of maximum diameter and minimum split help to improve the performance?

To answer the questions Q1 and Q2, we run both models, without giving any criterion or user-constraints, to nd all partitions for a clustering task with 2 clusters. Both models CP1 and CP2 are implemented on Gecode solver version 4.2.1. The experiments are performed on a 3.4GHz Intel Core i5 processor running Ubuntu 12.04. They both use the same branching rule: in CP1, variables are chosen in order I 1 , . . . , I k then G 1 , . . . , G n whereas in CP2, variables are chosen in order G 1 , . . . , G n . For the choice of values for each variable, the values are enumerated in the ascending order. The choice of the dataset is not important as the dissimilarities between points are not used in partition constraints.

The number of partitions depends only on the number of points n and the number of clusters k.

Theorem 23 Given a number of points n and a clustering task that divides these points into 2 clusters, the number of dierent partitions, when ignoring the index of each cluster in a partition, is given by a function f (n), which is dened by the recursive relation:

f (n) = 1 if n = 2, 2 × f (n -1) + 1 if n > 2.
P = (C 1 = {o 1 }, C 2 = {o 2 })
In case n > 2, for each partition P = (C 1 , C 2 ) of n-1 points {o 1 , . . . , o n-1 } into 2 clusters, there are two corresponding partitions with n points {o 1 , . . . , o n }:

P 1 = (C 1 ∪ {o n }, C 2 ), P 2 = (C 1 , C 2 ∪ {o n }).
Moreover, with n points, we have one more partition: ({o 1 , . . . , o n-1 }, {o n }). In total, there are 2 × f (n -1) + 1 partitions of n points into 2 clusters.

If there are no symmetric solutions, the total number of solutions given by our models must be equal to the number of partitions given by the function f (). Table 6.1 shows the results of our experiments. We present in the rst column the number of points n and in the second column the value of f (n). The third and fourth columns present the number of solutions and the total run time taken by CP1. The next two columns express the results corresponding to CP2. The last column expresses the speed up from CP1 to CP2, given in percentage.

It is clear that both our models nd the same number of solutions as the value of f (n). It shows that, in practice, there is no symmetry solutions found by both models. The model CP2 is 60% faster than CP1. In all the tests, the number of nodes in the search tree is identical in both models. It shows that the modeling of variables and partition constraints in CP2 is more ecient than CP1.

To respond the question Q3, we take the same datasets and the same experimentation as presented in previous chapter in Table 5.1.

For the criterion of maximal diameter We consider the clustering task without user-constraints, and the number of clusters k is xed (for CP2:

k min = k max = k).
Similarly to the previous chapter, in the experiments, the time-out is set to 1 hour and the Euclidean distance is used to compute the dissimilarity between objects. The number of clusters k is set to the real number of classes given in Table 5.1. Table 6.2 shows the results of our experiments. We present in the rst column the datasets. The second and • The modeling of variables and partition constraints: the improvement is clear as explained in question Q 2 .

• The heuristic of branching: in CP2, for selecting the value for branching, we consider all assigned variables, whereas in CP1, only the representative point is measured. This improvement is signicant, in some datasets (for example: User Knowledge, Vehicle, Yeast, Image Segmentation) there are much less number of nodes in the search tree generated by CP2 than CP1. As a result, the total run time of CP2 in those datasets is much faster than CP1. Figure 6.3 compares the diameter of the rst solution found by both CP1 and CP2. We compute the ratio of the diameter found in the rst solution to the optimal Diameter. With a good heuristic on branching, the ratio should be close to 1. CP2 has a better ratio in most datasets, except for the dataset Synthetic Control. In many tests, the ratio in CP2 is less than 1.1, which means that the diameter of the rst solution is very close to the optimal diameter. It explains the boost on performance with the second model. For the criterion of split We redo the experiments of split done in the previous chapter, with the second model (CP2). We consider the clustering task with the criterion of split and with an user-constraint: the diameter of each cluster must not exceed 1.5D opt , where D opt is the optimal diameter. In addition, we experiment the second model in two cases: the number of classes k is xed to the actual number of classes, and the number of classes k is bound from 2 to the actual number of classes. Table 6.3 presents the results. Similar to the criterion of diameter, the performance of CP2 is much better than that of CP1. There is a slight increase in time taken in the second case, as the model has to consider more partitions.

Experiments on user-constraints

We evaluate the eectiveness of our constrained clustering model (CP2), in term of the quality of the optimal partition with user-constraints and a given criterion. For measuring a partition, we consider the Adjusted Rand Index (ARI) [START_REF] Hubert | [END_REF]. The ARI measures the similarity between two partitions, in this case, the true partition P of the dataset and the partition P found by our model. The ARI is dened as: where a is the number of pairs of points that are in the same cluster in P and in P , b is the number of pairs of points that are in dierent clusters in P and in P , c is the number of pairs of points that are in the same cluster in P , but in dierent clusters in P and d is the number of pairs of points that are in dierent clusters in P , but in the same cluster in P . We consider the instance-level constraints: must-link and cannot-link constraints. The constraints are generated following the method in [START_REF] Wagsta | [END_REF]: two points are chosen randomly from the dataset, if they are in the same cluster in P , a must-link constraint is generated, otherwise a cannot-link constraint is generated. For the criterion, we consider the partition given by: a. the rst solution with the criterion of Diameter.

ARI = 2(ab -cd) (a + d)(d + b) + (a + c)(c + b)
b. the optimal solution with the criterion of Diameter. c. the rst solution with the criterion of WCSD. d. the optimal solution with the criterion of Split.

The experiments are on the two datasets Iris and Wine. Surprisingly, the criterion of split cannot solve these datasets with pairwise constraints. The underlying reason is that in these datasets, there are two clusters that are very close. As a result, the minimum split in the real partition is close to 0. The propagation of the constraints for the criterion of split is poor and the search space is too large. The solver has to consider many partitions and cannot end the search in one hour. We report in Figure 6.4 and 6.5 the average ARI over 100 constraint sets for 0 (unconstrained), 30, 60, . . . , 300 constraints randomly generated from the data. In both datasets, the ARI is improved as the number of constraints increases. For the dataset Wine, the ARI of partitions found without userconstraint is typically worse. However, the ARI improves when there are more than 200 constraints (200 is equivalent to 1.27% of the number of pairs of points in this dataset). We evaluate the inuence of the bounds on the number of clusters k ∈ [k min , k max ] on the performance with the diameter criterion. In the rst experiment, we set k max = 10 and we increment k min from 2 to 10. The diameter criterion favours the high number of clusters such that the higher the number of clusters, the smaller the value of the optimal diameter. Without user-constraints, the optimal solution with the criterion of diameter always has a number of clusters equal to k max . The result is given in Figure 6.6 where we report the number of nodes in the search tree, for each dataset, when we vary the bound k min . For all datasets, the number of nodes in the search tree is constant when k min is changed. By changing the bound, the performance of our model remains the same. It shows that the propagation of the constraint of diameter is eective. After nding and proving the optimal solution with k max clusters, the solver can conclude that there does not exist a better solution, with less than k max clusters.

In the second experiment, we x k min = 2 and we increment k max from 2 to 10. The result is given in Figure 6.7 where we report the number of nodes in the search tree, for each dataset, when we increment the bound k max . In general, when k max increases, there are more partitions to be considered. However, we see an interesting trend in Figure 6.7. In many datasets, the number of nodes in the search tree does not increase as k max increases. Indeed, since k max is higher, the optimal maximum diameter is smaller, and propagation of the diameter constraint is more eective. It explains that in some cases, it takes less computation time when k max increases.

Summary

In this chapter, we proposed a second model in Constraint Programming for solving the problem of Constrained Clustering. Its goal is to maximally improve the performance as well as to give more exibility to users . Our rst model presented in Chapter 5 is based on a two-level representation: a set of variables for the assignment of a representative for each class and a set of variables for the assignment of a representative for each point. This choice of variables requires that the number of classes k is xed beforehand, since each representative is modelled by a CP variable. The second model is based only on a set of variables for the assignment of a number of class to each point. In this model, the number of classes k is only bounded by k min and k max given by the user. We designed two new global constraints to replace a large number of reied constraints for the criteria of diameter and split. The dierence between two models is signicant. We demonstrated that while the second model is more modular and oers more exibility, it is more ecient than the rst model. In real world, many problems are concerned with dierent conicting criteria. For instance, when we consider to buy a new car, we wish to get a car that has high speed, high durability, low petrol consumption, but with low price. The four criteria are clearly in conict and it is impossible to nd such a car that is the best in each criterion whereas nding the best car with any of the four criteria alone is easy. This is an example of a multi-objective optimization problem. In mathematics, multi-objective optimization is the problem of optimization where two or more objective functions are considered simultaneously. A Multi-objective Optimization Problem (MOP) is dened as follows [?]:

M inimize f (x) = [f 1 (x), . . . , f m (x)] T subject to x ∈ X
where m is the number of criteria. The vector x ∈ R n is composed of n decision variables. Each vector x, also called point x, presents a solution in the multi-objective optimization problem. X is the decision space, also called the feasible design space, which describes the feasible solutions that can be made. The criterion space Z is dened as the set {f (x)|x ∈ X }. The vector function f : R n → R m is formed by m objective functions f i : R n → R. In the case where some objective functions are to be maximized, it is equivalent to minimize its negative.

Figure 7.1 expresses an example of the decision space X and the criterion space Z. Note that a point in the criterion space may correspond to more than one point in the decision space.

Typically, the objective functions conict each others and there is no single optimal solution that optimizes all objective functions. For evaluating solutions of MOPs, it is necessary to dene how to compare dierent objective vectors f (x) for dierent solutions x ∈ X . A common way to dene it is using the Pareto dominance relation [START_REF] Pareto | Cours d'economie politique[END_REF]].

Denition 24 Given two vectors of function values f (x 1 ) and f (x 2 ) ∈ Z, w.r.t two solutions x 1 and x 2 . We say the vector f (x 1 ) Pareto-dominates f (x 2 ), denoted by f (x 1 ) ≺ f (x 2 ), if and only if the solution x 1 is no worse than x 2 in all objectives and the solution x 1 is strictly better than x 2 in at least one objective.

f (x 1 ) ≺ f (x 2 ) i f i (x 1 ) ≤ f i (x 2 ) ∀i ∈ 1, . . . , m, ∃j ∈ 1, . . . , m f j (x 1 ) < f j (x 2 ).
For simplication, we also say that solution x 1 dominates solution x 2 , denoted by

x 1 ≺ x 2 , if f (x 1 ) ≺ f (x 2 ).
From the denition of the Pareto dominance relation, it is obvious that the "best" solutions of a MOP should not be dominated by any other solutions. The solutions that are not dominated by any other solutions are called Pareto optimal solution, dened as follows:

Denition 25 A feasible solution x * ∈ X is called Pareto optimal, if there is no other

x ∈ X such that f (x) ≺ f (x * ).
In practice, a common way to solve a MOP is nding the Pareto optimal set, or the Pareto Front set.

Denition 26 The Pareto optimal set of solutions P * is dened as follows:

P * = {x ∈ X | y ∈ X : f (y) ≺ f (x)}
Denition 27 The Pareto Front set PF * of a MOP is the projection of its Pareto optimal set in the criterion space.

PF * = {f (x) | x ∈ P * }
Example 27 Let us consider a MOP with two objective functions f 1 and f 2 to be minimized. Figure 7.2 illustrates an example of Pareto optimal set and Pareto front set. In this example, solution x 1 does not dominate any other solutions, as f 2 (x 1 ) is the maximized value. However, the solution x 1 is Pareto optimal, as it is not dominated by any other solutions. The solution x 3 is not Pareto optimal. The reason is that f 1 (x 3 ) > f 1 (x 2 ) and f 2 (x 3 ) > f 2 (x 2 ), so solution x 3 is dominated by solution x 2 .

Bi-criterion Constrained Clustering

Clustering with the criterion of minimizing the maximal diameter aims at nding homogeneous clusters, but it often suers from the dissection eect [Cormack 1971], i.e. quite similar objects may be classied in dierent clusters, in order to keep the diameter small. On the other hand, clustering with the criterion of maximizing the minimal split, which aims at nding well separated clusters, often suers from the chain eect [Johnson 1967], i.e. a chain of close objects may group very dierent objects in the same cluster. The popular WCSS criterion minimizes the sum of the squared distances between points and the center of their cluster. With this criterion sometimes objects which should be in a large group may be classied in dierent clusters in order to keep this sum small. An ideal partition of homogeneous and well-separated clusters should have a minimum diameter and a maximum split. Unfortunately, such a partition in general does not exist, since the two criteria are often conicting. [Delattre & Hansen 1980] suggests to consider the bi-criterion optimization problem of maximizing the minimal split between clusters and minimizing the maximal cluster diameter. This bi-criterion is natural and captures both the homogeneity and the separation requirements for a good clustering. [Delattre & Hansen 1980] proposes an algorithm for searching the complete Pareto front set, w.r.t these two criteria. In this thesis, we consider the same bi-criterion problem, but with the user-constraints. As presented in section 7.1, a Pareto optimal solution is a solution for which it is impossible to improve the value of one criterion without degrading the value of the other criterion. Considering the two criteria of diameter and split, given two partitions ∆ and ∆, we say ∆ dominates ∆, denoted by ∆ ≺ ∆, if and only if:

D(∆ ) ≤ D(∆) and S(∆ ) > S(∆) or D(∆ ) < D(∆) and S(∆ ) ≥ S(∆).
A partition ∆ is Pareto optimal if and only if there is no other partition ∆ that dominates ∆. If the user species a function on the criteria to optimize, for example max(S/D) or min[αD -(1 -α)S] with 0 ≤ α ≤ 1, the optimal solution will be among the Pareto optima [START_REF] Wang | [END_REF]. For instance, when the number of classes is set to 3, the Pareto optima of the example in Figure 7.3 are given in Figure 7.4. If we minimize the ratio S/D, the optimal solution is the solution 5, which seems the most reasonable. The user can specify polynomial algorithm is proposed in [START_REF] Wang | [END_REF][START_REF] Wang | [END_REF]. However, to the best of our knowledge, there is no algorithm dealing with this bi-criterion, while supporting any kind of user constraints.

CP Model for Bi-criterion Clustering with userconstraints

Our CP model represents a general and declarative framework for constrained clustering, where a user can choose one among dierent optimization criteria and can integrate different kinds of user constraints. This exibility oers dierent possibilities of usage. We present in this section the use of our model to deal with bi-criterion constrained clustering tasks. We are given a constrained clustering task with a set of user constraints U, which is possibly empty. We aim at computing the Pareto front set for this constrained clustering task with the bi-criterion (max S, min D). One approach to achieve this, is given in Algorithm 15; it is comparable to the -constraint approach as presented in [T'kindt & Billaut 2006]. In this algorithm, optimizations with a single criterion are iterated, each time with a condition on the value of the other criterion. The function Maximize_Split(U) or Minimize_Diameter(U) means the use of our model with the criterion of maximizing the split or minimizing the diameter, respectively, and with the set of user constraints U. It returns an optimal solution which satises all the constraints in U, if there exists one, or N U LL otherwise. The idea of the algorithm is presented in Figure 7.5. The rst Pareto solution found is the solution with the best diameter whereas the last Pareto solution found is the solution with the best split. We prove that this algorithm nds the complete Pareto front set.

A ← ∅; i ← 1; ∆ D i ← Minimize_Diameter(U); while ∆ D i = N U LL do ∆ S i ← Maximize_Split(U ∪ {D ≤ D(∆ D i )}) ; A ← A ∪ {∆ S i }; i ← i + 1; ∆ D i ← Minimize_Diameter(U ∪ {S > S(∆ S i-1 )}) ;
. if ∆ D i = N U LL then ∆ S i = N U LL,
2. for all 1 < i ≤ m, S(∆ S i ) > S(∆ S i-1 ),

for all

1 < i ≤ m, D(∆ D i ) > D(∆ D i-1 ), 4. for all 1 ≤ i ≤ m, D(∆ S i ) = D(∆ D i ),
5. ∃∆ such that D(∆) < D(∆ D 1 ),

for all

1 ≤ i < m, ∃∆ such that S(∆) ≥ S(∆ S i ) and D(∆) < D(∆ D i+1 ). Proof 1. If ∆ D i = N U LL, since ∆ D i satises the constraint D ≤ D(∆ D i ), it is a candidate for ∆ S i .
2. Since among all the partitions satisfying D ≤ D(∆ D i ), ∆ S i is the one which maximizes the split (line 5), we have S(∆ S i ) ≥ S(∆ D i ). Since ∆ D i is among the partitions satisfying S > S(∆ S i-1 ) (line 8), we have S(∆ D i ) > S(∆ S i-1 ). Therefore S(∆ S i ) > S(∆ S i-1 ).

3. The set of the partitions satisfying S > S(∆ S i ) is a strict subset of the set of the partitions satisfying S > S(∆ S i-1 ), since S(∆ S i ) > S(∆ S i-1 ). Since ∆ D i and ∆ D i-1 are the partitions which minimize the diameter among all the ones in these two respective sets, we have

D(∆ D i ) ≥ D(∆ D i-1 ). If D(∆ D i ) = D(∆ D i-1 ), then S(∆ S i ) = S(∆ S i+1 ) (line 5) which contradicts S(∆ S i ) > S(∆ S i-1 ). Therefore D(∆ D i ) > D(∆ D i-1 ).
S i ) ≤ D(∆ D i ). Since S(∆ S i ) > S(∆ S i-1
) and ∆ D i is a partition which minimizes the diameter among all the ones satisfying the condition S > S(∆ S i-1 ), we have

D(∆ S i ) ≥ D(∆ D i ). Therefore D(∆ S i ) = D(∆ D i ),
5. Since ∆ D 1 is a partition which minimizes the diameter under no condition (line 3), there exists no partition ∆ with D(∆) < D(∆ D 1 ). ). The rst case is impossible since ∆ S i is a partition which maximizes the split among all those which satisfy the condition D ≤ D(∆ D i ). For the second case, if i = 1 then ∆ does not exist according to the point 5 of Proposition 28. If i > 1, since S(∆ S i ) > S(∆ S i-1 ), the partition ∆ must satisfy D(∆) < D(∆ D i ) and S(∆) > S(∆ S i-1 ). This is impossible according to the point 6 of Proposition 28.

2. Let us assume that there exists a Pareto optimal solution such that (D(∆), S(∆)) ∈ B. According to the points 5 and 6 of Proposition 28, ∆ to be not dominated must have S(∆) > S(∆ S m ). But Algorithme 15 terminates after ∆ S m , i.e. there exists no partition with the split greater than S(∆ S m ). Therefore, a Pareto optimal solution ∆ ∈ A does not exist.

Multi-objective optimization in Constraint Programming in only one phase of search is proposed in [Gavanelli 2002]. The idea is to realize a global constraint P areto(Obj 1 , . . . , Obj m , A), which keeps a set of non dominated solutions so far computed A and which operates on the variables representing the objective functions Obj i . This constraint reduces the domain of a variable Obj i if the domains of the other variables enter into the dominated zone of a solution in A. A detailed description of this constraint Pareto can be introduced in our model. Nevertheless this approach for the moment is still much less ecient than Algorithm 15 and therefore needs deeper studies, for instance of the search strategy, in order to improve the eciency.

Experimentation

Performance test

For the bi-criterion split-diameter, we compare our second model presented in Chapter 6 (denoted by CP2) with the bi-criterion clustering algorithm based on graph coloring [Delattre & Hansen 1980] (denoted bGC). The program is not available, we have coded it in C++. To our knowledge, this is the only exact algorithm for general values of k.

In the experiments, datasets presented in Table 5.1 are used and the time-out is set to 1 hour. The number of clusters k varies between 2 and the real number of classes. Table 7.1 shows the result of the experiments. The second column (#Sol) gives the number of Pareto optimal solutions found and the next columns give the run time in seconds of each approach. Note that both approaches nd a complete Pareto front set, so the number of Pareto optimal solutions found is the same as it is the cardinality of the complete Pareto front set. However, the partitions found by the two approaches might be dierent as there exists dierent partitions with the same values of diameter and split. It can be seen that our model is the most ecient in most cases. It takes advantage of the ecient constraint propagation to reduce the search space. As in the case of GC, the algorithm bGC is limited to datasets with less than 500 points.

Bi-criterion clustering with user-constraints

We consider the instance-level constraints: must-link and cannot-link constraints. The constraints were generated with the same method as in the previous Chapter: two points are chosen randomly from the dataset, if they are in the same cluster, a must-link constraint is generated, otherwise a cannot-link constraint. We generate a set of 80 instancelevel constraints from the dataset Iris and use our second model for the task of bi-criterion constrained clustering. The rst test is without user-constraints, the second one is with the rst 20 instance-level constraints, the third one is with the rst 40 instance-level constraints and so on. Figure 7.6 presents the results of the Pareto front in six cases, where there are from 0 to 80 instance-level constraints. As there are more and more userconstraints, it is clear that there are less feasible solutions and as a result, the criterion space is reduced. Because we generated user-constraints from real labels, it is obvious that the values (diameter,split) of the real partition of the dataset lies in the smallest criterion space (the region dened by red lines in Figure 7.6). There are signicantly dierence of the location of Pareto front between dierent cases in Figure 7.6. It shows that without user-constraints, even there are many points in the Pareto front, they may all lie far from the real partition in the criterion space. For that reason, it could be useful to enable user-constraints for the task of bi-criterion clustering. In this chapter we explore the problem of constrained clustering with the multi-objective optimization. Typically, for solving a multi-objective optimization problem, a common way is to nd the set of Pareto optimal solutions. We proposed an ecient algorithm, based on the -constraint approach [T'kindt & Billaut 2006], using our CP model for nding the 7.5. Summary 107 complete Pareto front set, w.r.t the two criteria of split and diameter. Experiments on various UCI benchmark datasets show that our model is more eective than the bi-criterion clustering algorithm based on graph coloring [Delattre & Hansen 1980]. More than that, to the best of our knowledge, we are the rst to integrate user-constraints to the problem of bi-criterion clustering. A notre connaissance, notre approche est la première pour les problèmes de classication non supervisée sous contraintes utilisateur, qui intègre le critère de diamètre ou de marge ou de la somme des dissimilarités intra-cluster et du bi-critère marge-diamètre, en présence de contraintes utilisateur.

Contributions of this Dissertation

Cluster analysis is an important task in Data Mining with hundreds of dierent approaches in the literature. Since the last decade, cluster analysis has been extended to constrained clustering, or semi-supervised clustering, so as to integrate previous knowledge on data to clustering algorithms. In this dissertation we dene two models based on Constraint Programming (CP) to address constrained clustering tasks. Relying on CP has two main advantages: the declarativity, which enables to easily add new constraints and the ability to nd an optimal solution satisfying all the constraints (when there exists one). However, modelling a problem in CP in an ecient and eective way requires a lot of skills and eorts.

The rst model is based on a two-level representation: a set of variables for the assignment of a representative for each class and a set of variables for the assignment of 110 Chapter 8. Conclusion and Future Work a representative for each point. In literature, most constrained clustering algorithms are limited to few kinds of user-constraints. Our model, based on CP, is more general in the sense that it supports instance-level constraints, specifying requirements on pairs of objects, as well as dierent cluster-level constraints, specifying requirements on clusters. It also allows user to choose among dierent optimization criteria. In order to improve the eciency, dierent aspects have been studied in the dissertation. We dene a ltering algorithm for the criterion of within cluster sum of dissimilarities in our model. Experiments on various UCI benchmark datasets show that our model is competitive with other exact approaches.

The second model is based only on a set of variables for the assignment of a cluster, represented by its number, to each point. The model is more exible as it does not require to x a priori the number of clusters k, but only bounds on k are required. The second model is lighter in terms of the number of variables, as well as the number of constraints. We dene two new global constraints for the criteria of split and diameter in order to replace a large number of reied constraints. Experiments show that the second model is more eective and can handle larger datasets.

Another important contribution is to apply the model to the task of constrained bicriterion clustering. We show that our model can easily be embedded into a more general process and we illustrate this on the problem of nding the Pareto front of a bi-criterion optimization process. To the best of our knowledge, we are the rst to:

• Find a global optimum for the problem of Constrained Clustering with the criteria of diameter, split and within cluster sum of dissimilarities.

• Solve the bi-criterion (split-diameter) clustering with user constraints.

Directions for Future Work

Directions for future work include the following:

• Firstly, we would like to study dierent heuristics for the search strategy in order to improve the performance of our models. The goal of the search strategy is not only to quickly nd an optimal solution, but to quickly prove the optimality. During the search process, the choices on variables and values for branching are important, as the correct choices may help the constraint propagation to reduce eciently the search space. A good search strategy may not only depend on the criterion, but also depend on the type of datasets. As a consequence, it is worth to take time and eort to study and to experiment with dierent techniques.

• The goal of our models is to nd an optimal solution for a given criterion. In practice, it is not always true that the optimal solution corresponds to a good partition whereas nding and proving an optimal solution requires exhaustive search. We would like to study and apply heuristics search techniques, to balance between the computation time and the quality of the solution. For instance, heuristics based on Large Neighbourhood Search [Shaw 1998] have recently shown tremendous results 8.2. Directions for Future Work 111 in solving various problems in CP. We would like to study and apply Large Neighbourhood Search for solving large-scale datasets, in which it is impossible to nd an optimal solution.

• Relating to the multi-objective clustering, our approach presented in Chapter 7 requires to restart the model multiple times in order to get the complete Pareto front. Multi-objective optimization in Constraint Programming in only one phase of search is proposed in [Gavanelli 2002]. The idea is to realize a global Pareto constraint, which keeps a set of non dominated solutions and operates on the variables representing the objective functions. A detailed description of this constraint as well as an extension with Large Neighbourhood Search is proposed in [?]. This Pareto constraint can be introduced in our model. Nevertheless this approach for the moment is still limited and much less ecient than the approach presented in Chapter 7. It therefore needs deeper studies, for instance of the ltering algorithm or of the search strategy, in order to improve the eciency.

• Our models nd an optimal solution that satises all given user-constraints. However, in practice, there may exist noise in the set of user-constraints, or the problem is over-constrained, i.e. there is no solution satisfying all constraints. We would like to integrate soft constraints to our model, i.e. constraints that are preferred but not required to be satised. Soft constraints allow users to set a threshold on the number of constraints that must be satised, while there are some constraints unsatised.

One solution is to study other modelings for the problem, for instance to design new dedicated ltering algorithms for soft constraints. Another possible technique that we would like to study is applying a schema that allows hard constraints for representing soft constraints [Régin 2011].

• Our models support the criteria of diameter, split and within cluster sum of dissimilarities. We would like to study and integrate other well-known criteria to the model, for instance, the within-cluster sum of squares. Adding a new criterion in our model is simple, by adding new constraints to express the function objective. There are dierent ways to model a given criterion in CP. The challenge is how to model a criterion with high performance. The propagation of criterion constraints is therefore highly important. We need ltering algorithms that have low computational complexity but are able to compute good bounds on the objective function and to lter as much as possible inconsistent values in the domain of variables. We would like to study and to apply techniques in linear programming, for instance column generation techniques, for designing eective ltering algorithms for other criteria.

• Finally, we would like to extend our research to address other clustering tasks, for instance fuzzy clustering where points can belong to more than one cluster. Mots clés : classication non supervisée, contraintes utilisateur, Programmation par Contraintes, clustering bi-critère.

Constrained Clustering by Constraint Programming

Abstract: Cluster analysis is an important task in Data Mining with hundreds of dierent approaches in the literature. Since the last decade, the cluster analysis has been extended to constrained clustering, also called semi-supervised clustering, so as to integrate previous knowledge on data to clustering algorithms. In this dissertation, we explore Constraint Programming (CP) for solving the task of constrained clustering. The main principles in CP are: (1) users specify declaratively the problem in a Constraint Satisfaction Problem; (2) solvers search for solutions by constraint propagation and search. Relying on CP has two main advantages: the declarativity, which enables to easily add new constraints and the ability to nd an optimal solution satisfying all the constraints (when there exists one). We propose two models based on CP to address constrained clustering tasks. The models are exible and general and supports instance-level constraints and dierent cluster-level constraints. It also allows the users to choose among dierent optimization criteria. In order to improve the eciency, dierent aspects have been studied in the dissertation. Experiments on various classical datasets show that our models are competitive with other exact approaches. We show that our models can easily be embedded in a more general process and we illustrate this on the problem of nding the Pareto front of a bi-criterion optimization process.

  Contributions . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 List of publications . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Dissertation Organization . . . . . . . . . . . . . . . . . . . . . . . 4

•

  Thi-Bich-Hanh Dao, Khanh-Chuong Duong and Christel Vrain. A Declarative Framework for Constrained Clustering. In Proceedings of the European Conference on Machine Learning and Principles and Practice of Knowledge Discovery in Databases, pages 419434, 2013 • Thi-Bich-Hanh Dao, Khanh-Chuong Duong and Christel Vrain. A Filtering Algorithm for Constrained Clustering with Within-Cluster Sum of Dissimilarities Criterion. In Proceedings of the 25th International Conference on Tools with Articial Intelligence, pages 10601067, 2013 National Conference Papers • Thi-Bich-Hanh Dao, Khanh-Chuong Duong and Christel Vrain. Une approche en PPC pour la classication non supervisée. In 13e Conférence Francophone sur l'Extraction et la Gestion des Connaissances EGC, 2013 • Thi-Bich-Hanh Dao, Khanh-Chuong Duong and Christel Vrain. Un modèle général pour la classication non supervisée sous contraintes utilisateur. In Neuvième Journées Francophones de Programmation par Contraintes, 2013 • Thi-Bich-Hanh Dao, Khanh-Chuong Duong and Christel Vrain. Classication non supervisée mono et bi-objectif sous contraintes utilisateur par la programmation par contraintes. In Dixième Journées Francophones de Programmation par Contraintes, 2014.

  Figure 2.1: graph value (left) and a maximum matching (right) of Example 12

Figure 2

 2 Figure 2.2: residual graph G M (left) and marked edges (right) in Example 13

Figure 2

 2 Figure 2.3: A search tree for a CSP

Figure 2 . 4 :

 24 Figure 2.4: Search tree with the strategy: the variable with the smallest domain is selected for instantiation

Figure 2 . 7 :

 27 Figure 2.7: Solution symmetries on 8-queens problem

  k-Means algorithm Input: dataset O = {o 1 , . . . , o n }, a number of clusters k Output: a partition P with k clusters C 1 , C 2 , . . . , C k 1 Generate randomly k objects as initial cluster centers: m 1 , m 2 , . . . , m k i is assigned to the closest centroid: 5 Assign o i to C c such that d(o i , m c ) = min j∈[1,k] d(o i , m j ) change of cluster centroids; 12 Return partition

  k-Medoids algorithm for partitioning Input: dataset O = {o 1 , . . . , o n }, a number of clusters k Output: a partition P with k clusters C 1 , C 2 , . . . , C k 1 Randomly generate k objects as initial cluster medoids x 1 , x 2 , . . . x k i is assigned to the closest medoid: 5 Assign o i to C c such that d(o i , x c ) = min j∈[1,k] d(o i , x j ) S when swapping the object o i with a representative x 9 if S < 0 then 10 o i replace x as a new medoid 11 end 12 until no change of cluster medoids; three dimensional space. Algorithm 4: FPF algorithm for partitioning Input: dataset O = {o 1 , . . . , o n }, a number of clusters k Output: a partition P with k clusters C 1 , C 2 , . . . , C k

6

  Compute the number of colors N color needed to colour the graph G 7 if N color ≤ k then 8 d max ← the next largest dissimilarity 9 end 10 until number of colors needed > number of clusters k;

Algorithm 7 :

 7 Branch-and-bound algorithm for clustering with Depth-First-Search 1 A heuristic is used to obtain a good solution 2 The upper bound is set to the criterion value of the solution 3 number of assigned objects p ← 0 c ∈ [1, k] do 13 Assign object o p to cluster C c 14 Evaluate the partial solution with p objects 15 if the partial solution may be a part of the optimal solution

  clusters with a branch-and-bound algorithm. The search strategy assigns objects into clusters in the order: o n-k-i+1 , . . . , o n . Considering a partial solution with p assigned objects, it is clear that the assigned objects are {o n-k-i+1 , . . . , o n-k-i+p } and the unassigned objects are {o n-k-i+p+1 , . . . , o n }. The optimal wcsd of objects {o n-k-i+p+1 , . . . , o n } into k clusters is obviously already computed in one of the previous steps, therefore the third part can be obtained optimally and the lower bound of the partial solution is evaluated better. Figure3.1 illustrates the steps of RBBA algorithm for solving the wcsd clustering with 7 objects o 1 , . . . , o 7 into 2 clusters.

Figure 3 . 1 :

 31 Figure 3.1: Illustration of steps of repetitive branch-and-bound algorithm

  o i do 4 Assign object o i to the closest cluster C c such that:

  (4.1) to a smaller set of columns T ⊆ T . The dual of the master problem, a LP relaxation, can be solved by three dierent strategies[START_REF] Du | [END_REF] and the best strategy is the ACCPM interior point method of[Gon et al. 1992]. Suppose that solving the restricted master problem returns values for σ and λ i . The objective of the auxiliary problem is to nd a 50 Chapter 4. Declarative Data Mining column t with negative reduced cost. The column t with the smallest reduced cost is the solution of: auxiliary problem,[START_REF] Du | [END_REF] considers the problem as a hyperbolic program in 0-1 variables, and solves it by a branch-and-bound algorithm for unconstrained quadratic 0-1 optimization. The column generation framework iteratively solves the LP relaxation of the master problem and the auxiliary problem to nd the optimal solution of the LP relaxation. This solution is not necessarily integer and branching is added to obtain an optimal partition. The column generation framework[START_REF] Du | [END_REF] solved for the rst time the optimal sum of squares partition on datasets with 100-200 objects, as for example the dataset Iris.[Aloise etal. 2012] improved the performance of the framework by introducing a new method based on geometric arguments to solve the auxiliary problem. It is claimed to solve the optimal partition with datasets over 2300 objects.[START_REF] Babaki | [END_REF] introduced for the rst time user-constraints to the framework.

Example 19

 19 Given a dataset of 7 points O = {o 1 , . . . , o 7 } and a clustering task dividing the points into 2 clusters. We dene I = {I 1 , I 2 } and G = {G 1 , . . . , G 7 } where the initial domains of all variables are {1, . . . , 7}. Let us consider a partition of 2 clusters, the rst one composed of o 1 , o 2 , o 5 and o 6 and the second one composed of the remaining points. The points are denoted by their integer indices (o 1 is denoted by 1, o 2 by 2 and so on).

Figure 5 .

 5 1 presents this example. The left image shows the dataset and the domain of variables at the beginning whereas the right image presents a partition and the corresponding assignment of variables.

Figure 5 . 1 :

 51 Figure 5.1: Example of assignment values to variables in a solution

Figure 5 . 2 :

 52 Figure 5.2: Example of symmetric solutions

  It requires each cluster to have a number of objects inferior to a predened threshold β. Similarity, for each value c ∈ {1, . . . , k}, we post a AtMost constraint:AtM ost(G, I c , β).

Figure 5 . 3 :

 53 Figure 5.3: Example of ltering

  1 A group is like a clique and the number of pairwise connections is the number of edges in the clique.

Figure 5 . 4 :

 54 Figure 5.4: Example ofV.lb = V 1 + V 2 .lb + V 3 .lb

Figure 5 . 6 :

 56 Figure 5.6: Example of search strategy for the criterion of WCSD

Figure 5 . 7 :

 57 Figure 5.7: Before and After reordering points

  Figure 5.9: Datasets

Figure 5 .

 5 Figure 5.10: Maximal diameter optimization

Figure 6 .

 6 Figure 6.1: Assignment of values to variables

  It requires each cluster to have a number of objects inferior to a predened threshold β. For each value c ∈ {1, . . . , k max }, we post a constraint AtMost: AtM ost(G, c, β).

  noted by diameter(G, D, d), which exploits the dissimilarity measure d between any two points and which operates on the set of variables G and the variable D. The constraint diameter(G, D, d) ensures that D is the maximal diameter of the clusters formed by the variables G 1 , . . . , G n . This means, for every couple 1 ≤ i < j ≤ n, this constraint ensures the constraints in(6.1). By developing the constraint diameter(G, D, d), we maintain all these relations in one constraint, which makes the model much less heavy. The ltering algorithm is presented in Algorithm 13, where Dom(D) = [D.min, D.max).

Figure 6 . 3 :

 63 Figure 6.3: Comparison of CP1 and CP2 with the criterion of Diameter on 1st solution

Figure 6 . 5 :Figure 6 . 6 :

 6566 Figure 6.4: Dataset Iris with user-constraints

Figure 6 . 7 :

 67 Figure 6.7: Experiments on bound of k: k min = 2, k max = 2, . . . , 10

Contents

Figure 7 .

 7 Figure 7.1: Illustration of the decision and the criterion space

   gives an illustration of these eects. Image A shows three groups that are easily recognized. Image B shows the obtained solution with the diameter criterion when the

Figure 7 .

 7 Figure 7.2: Illustration of Pareto optimal set and Pareto front set

Figure 7 .

 7 Figure 7.3: Illustration of eects with dierent criteria

Figure 7 . 4 :

 74 Figure 7.4: Pareto optimal solutions

end 7 .Figure 7 . 5 :

 775 Figure 7.5: The solutions found by Algorithm 15

6 .

 6 Assume that there exists a solution ∆ such that S(∆) ≥ S(∆ S i ) and D(∆) < D(∆ D i+1 ). Since S(∆) ≥ S(∆ S i ) and ∆ D i+1 is a partition which minimizes the diameter among all those which satisfy the condition S > S(∆ S i ), we haveD(∆) ≥ D(∆ D i+1). This contradicts the fact that D(∆) < D(∆ D i+1 ).By Proposition 28, the position of the solutions found by Algorithm 15 is as presented in Figure 7.5. There exists no solution in the grey zone. A solution in the white zone is dominated by a solution ∆ S i ∈ A. Proposition 29 Given the set A = {∆ S 1 , . . . , ∆ S m } computed by Algorithm 15, the Pareto front set B = {(D(∆ S 1 ), S(∆ S 1 )), . . . , (D(∆ S m ), S(∆ S m ))} is complete, i.e.: 1. ∆ S i (1 ≤ i ≤ m) is a Pareto optimal solution, 2. ∃∆ Pareto optimal solution such that (D(∆), S(∆)) / ∈ B Proof Algorithm 15 terminates since according to Proposition 28, S(∆ S i ) > S(∆ S i-1 ) and these values are discrete and limited by the maximal dissimilarity of pairs of points. 1. We prove that for all ∈ [1, m], there exists no partition ∆ which dominates ∆ S i , i.e. D(∆) ≤ D(∆ S i ) et S(∆) > S(∆ S i ), or D(∆) < D(∆ S i ) et S(∆) ≥ S(∆ S i ). Since D(∆ S i ) = D(∆ D i ), the partition ∆ must satisfy D(∆) ≤ D(∆ D i ) et S(∆) > S(∆ S i ), or D(∆) < D(∆ D i ) et S(∆) ≥ S(∆ S i

  Figure 7.6: Bi-criterion constrained clustering with dataset Iris

  Dans ce chapitre, nous rappelons nos contributions et discutons quelques perspectives. Les deux premières contributions sont deux modèles basés sur la Programmation par Contraintes pour le clustering sous contraintes utilisateur. Les modèles sont généraux et exibles, ils permettent d'intégrer des contraintes d'instances must-link et cannot-link et diérents types de contraintes sur les clusters. Ils orent également à l'utilisateur le choix entre diérents critères d'optimisation. An d'améliorer l'ecacité, divers aspects sont étudiés. Les expérimentations sur des bases de données classiques et variées montrent qu'ils sont compétitifs par rapport aux approches exactes existantes. La dernière contribution de la thèse est de montrer que nos modèles peuvent être intégrés dans une procédure plus générale et nous l'illustrons par la recherche de la frontière de Pareto dans un problème de clustering bi-critères sous contraintes utilisateur.

Contents 8 . 1

 81 Contributions of this Dissertation . . . . . . . . . . . . . . . . . . 109 8.2 Directions for Future Work . . . . . . . . . . . . . . . . . . . . . . 110 This chapter concludes the dissertation and provides perspectives for future work.

  , Dao et al. 2013c, Dao et al. 2013a, Dao et al. 2013b].Chapter 6 presents our second model for constrained clustering. The model is more exible as it does not require to set a priori the number of clusters k, but only bounds on k are required. The second model is lighter in terms of the number of variables, as well as the number of constraints. The contributions of this Chapter are in[Dao et al. 2014a, Dao et al. 2014c, Dao et al. 2014b].

	1.3. Dissertation Organization	5

  Constraint Propagation . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Global constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.4 Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Search for COPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.5 Symmetries in CP . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

		.	15
	2.4.2	Variable and value ordering heuristics . . . . . . . . . . . . . . . . .	16
	2.4.3

Contents 2.1 General Concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 2.4.1 Backtracking search in CP . . . . . . . . . . . . . . . . . . . . . . .

  , IBM ILOG CP Optimizer [Laborie 2009], . . .

	Chapter 2. Constraint Programming

  ].

		Chapter 3
		Clustering
	Ce chapitre est consacré à l'état de l'art de la classication non supervisée (clustering) et
	du clustering sous contraintes, notamment les approches concernant nos recherches. Nous
	présentons les notions principales et diérents critères pour le problème de clustering. Ce
	chapitre introduit ensuite des méthodes fondées sur le partitionnement des données et
	des méthodes hiérarchiques. An de mieux modéliser la tâche d'apprentissage, mais aussi
	dans l'espoir de réduire la complexité, des contraintes dénies par l'utilisateur peuvent être
	ajoutées. Plusieurs types de contraintes peuvent être considérés; ces contraintes peuvent
	porter sur les clusters, comme par exemple sur leur diamètre ou sur leur taille, ou porter
	sur des paires d'objets qui doivent être ou ne doivent pas être dans une même classe.
	Nous présentons ces contraintes et les approches développées pour intégrer des contraintes
	utilisateur dans la tâche de clustering.
	Contents
	3.1 Problem denitions . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
	3.2 Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
	3.2.1	Criteria for a cluster . . . . . . . . . . . . . . . . . . . . . . . . . . .	25
	3.2.2	Criteria for a partition . . . . . . . . . . . . . . . . . . . . . . . . . .	27
	3.3 Partitioning methods . . . . . . . . . . . . . . . . . . . . . . . . . . 28
	3.3.1	Heuristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	28
	3.3.2	Exact approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	33
	3.4 Hierarchical clustering methods . . . . . . . . . . . . . . . . . . . . 36
	3.5 Constrained Clustering . . . . . . . . . . . . . . . . . . . . . . . . . 37
	3.5.1	Clustering algorithms for instance-level constraints . . . . . . . . . .	39
	3.5.2	Clustering algorithms for cluster-level constraints . . . . . . . . . . .	41
	3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
	Cluster analysis is an important task in Data Mining, which aims at partitioning a
	given set of objects into homogeneous and well-separated subsets called classes or clusters.

  chapitre est consacré à l'état de l'art des approches déclaratives en fouille de données, notamment les approches récentes concernant le domaine du clustering. Un bref aperçu des approches déclaratives est donné. Nous introduisons ensuite les travaux utilisant la programmation par contraintes [Khiari et al. 2010, Guns et al. 2013] pour résoudre le problème de recherche de k motifs fréquents. Ces cadres sont appliquées au problème de clustering conceptuel. Ce chapitre introduit ensuite diérents travaux utilisant la satisabilité des formules booléennes (SAT) [Davidson et al. 2010, Metivier et al. 2011, Berg & Jarvisalo 2013] et la programmation linéaire [Babaki et al. 2014] pour résoudre la tâche de clustering sous contraintes. Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 A large number of approaches have been developed for solving data mining problems.

		Chapter 4
		Declarative Data Mining
	4.2.1	A Constraint Language based on SAT for Declarative Pattern Mining	45
	4.2.2	A SAT framework for constrained clustering with 2 clusters . . . . .	46
	4.2.3	MaxSAT framework for Optimal Constrained Correlation Clustering	47

Ce Contents 4.1 Constraint Programming framework for k-pattern set mining . . 44 4.2 SAT based framework for clustering . . . . . . . . . . . . . . . . . 45 4.3 Column generation framework for Constrained Sum of Squares Clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.4

  Density constraint: The -constraint introduced in[Davidson & Ravi 2005b] requires for each point o i to have in its neighborhood of radius at least another point belonging to the same cluster:∀c ∈ [1, k], ∀o i ∈ C c , ∃o j ∈ C c , o j = o i and d ij ≤. This constraint tries to capture the notion of density, introduced in DBSCAN[Ester et al. 1996]. We have extended it by proposing a density-based constraint, stronger than the -constraint: it requires that for each point o i , its neighborhood of radius contains at least m points belonging to the same cluster as o i , i.e. ∀c ∈ [1, k],

• Minimal size α of clusters: It requires that each cluster has a number of points greater than a given threshold

α, i.e. ∀c ∈ [1, k], |C c | ≥ α.

• Maximal size β of clusters: It requires each cluster to have a number of objects inferior to a predened threshold

β, i.e. ∀c ∈ [1, k], |C c | ≤ β.

• Split δ-constraint: It expresses that the split between any two clusters must be at least

δ, i.e. ∀c ∈ [1, k], ∀c = c, ∀o i ∈ C c , ∀o j ∈ C c , d ij ≥ δ.

• Diameter constraint γ: It expresses that the diameter of each cluster must be at most γ, i.e. ∀c ∈ [1, k], ∀o i , o j ∈ C c , d ij ≤ γ.

•

Table

	5.7. Experiments				75
	Dataset	# Objects # Attributes # Classes
	Iris	150		4	3
	Wine	178	13	3
	Glass	214		9	7
	Ionosphere	351	34	2
	User Knowledge	403		5	4
	WDBC	569	30	2
	Synthetic Control	600	60	6
	Vehicle	846	18	4
	Yeast	1484		8	10
	Multiple Features	2000	10	6
	Image Segmentation	2000	19	7
	Waveform	5000	40	3
	Table 5.1: Properties of datasets
	Dataset	D opt	BaB	GC	CP
	Iris	2.58	1.4	1.8	< 0.1
	Wine	458.13	2	2.3	0.3
	Glass	4.97	8.1	42	0.9
	Ionosphere	8.6	-	0.6	0.5
	User Knowledge	1.17	-	3.7	84.5
	WDBC	2377.96	-	1.8	0.7
	Synthetic Control	109.36	-	-	56.1
	Vehicle	264.83	-	-	14.9
	Yeast	0.67	-	-2388.1
	Multi Features	12505.5	-	-	-
	Image Segmentation	436.4	-	-	649.2
	Waveform	15.6	-	-	-

  Chapter 5. An initial CP Model for Constrained Clustering

	Dataset	S opt Total time
	Iris	0.53	0.2
	Wine	53.33	0.3
	Glass	1.78	1.3
	Ionosphere	5.29	52.1
	User Knowledge Modeling	0.32	423.1
	WDBC	421.99	7.1
	Synthetic Control	43.59	31.3
	Vehicle	27.06	21.9
	Yeast

Table 5 .

 5

		4: split constraint and must-link constraints with dataset Iris
	split constraint WCSD Total time	# must-link	WCSD Total time
	no constraint	573.552	4174	no constraint	573.6	4174
	δ = 2% d max	573.552	1452	0.2%	602.6	1275.1
	δ = 4% d max	573.552	84.4	0.4%	602.6	35.6
	δ = 6% d max	573.552	0.3	0.6%	617	16.1
	δ = 8% d max	2169.21	0.1	0.8%	622.5	3.5
	δ = 10% d max	2412.43	<0.1	1%	622.5	1.6
	δ = 12% d max	2451.32	<0.1	100%	646	<0.1

  In this algorithm, if D.max has been changed and if D.max < d ij , the relation (6.1) will infer G i = G j , but we lter the domain of G i (or G j ) only if the other variable has been instantiated. Otherwise, if some variables in G have been instantiated, if G i = G j , with the relation (6.1) we infer D ≥ d ij so the lower bound D.min can be changed. Let us notice that as soon as the domain of one variable becomes empty, a failure case is detected by the solver. The worst case complexity is O(n 2 ).

Table 6 .

 6 1: Comparison of the two models without optimization criterionProof It is obvious that in the case there are two points {o 1 , o 2 }, there is only one possible partition P with 2 clusters, where each cluster contains one point:

	90	Chapter 6. A Modular CP Model for Constrained Clustering
			CP1		CP2
	n	f (n) #Solutions Time #Solutions Time Performance gain
	18	131071	131071	0.5	131071	0.3	167%
	19	262143	262143	1.1	262143	0.7	157%
	20	524287	524287	2.3	524287	1.4	164%
	21	1048575	1048575	4.8	1048575	3	160%
	22	2097151	2097151	10	2097151	6.1	164%
	23	4194303	4194303	20.8	4194303	12.8	163%
	24	8388607	8388607	42.8	8388607	26.1	164%
	25 16777215	16777215	88.8	16777215	53.7	165%

Table 6 .

 6 2: Performance of the two models with the criterion of minimizing the maximal diameter third column present the number of nodes in the search tree and the total run time taken by CP1. The two next columns express the results corresponding to CP2. The last column expresses the performance gain from CP1 to CP2, given in percentage. It is clear that the model CP2 outperforms CP1 on all datasets. All twelve datasets can be solved by CP2 within one minute. Although our two models are based on the Constraint Programming framework, there are three reasons that explain the signicant dierences in performance of the two models.

	6.6. Experiments

•

  The dedicated global constraint for the criterion of diameter: The improvement is clearer with larger datasets. For example with the dataset Yeast, the number of nodes in CP1 is 13 times more than CP2 but CP2 is 459 times faster than CP1. The reason is, in CP1, for modeling the criterion, the number of reied constraints is a square function of the number of points n. Although our dedicated global constraint

	92		Chapter 6. A Modular CP Model for Constrained Clustering
		1.6						1.59
	ratio to the optimal Diameter	1 1.2 1.4	1.17 1.04 1.02 1.03 1.02 1.09 1.33 1.03	1.12 1.09	1	1	1.09 1.15	1	1.36 1.1	1.28 1.01
			Iris Wine Glass Iono User WDBC Synt Vehi Yeast Image
				CP1 CP2

Table 6 .

 6 3: Comparison with the criterion of split

			Case 1		Case 2
		S opt	CP1	CP2	S opt	CP2
	Iris	0.53	0.2	< 0.1 0.53	< 0.1
	Wine	53.33	0.3	< 0.1 53.33	< 0.1
	Glass	1.78	1.3	0.1	1.78	0.7
	Ionosphere	5.29	52.1	1.2	5.29	1.2
	User Knowledge	0.32	423.1	2.8	0.32	2.9
	WDBC	421.99 7.1	0.6	421.99 0.6
	Synthetic Control	43.59	31.3	6.4	45.2	6.5
	Vehicle	27.06	21.9	3.7	27.06	3.7
	Yeast	0.15	-	106.3 0.15	106.3
	Multi Features	1107.1 -	82.2	1107.1 82.9
	Image Segmentation 228.7	-	20.4	-	20.6
	Waveform	-	-	-	-

  , nous montrons que nos modèles peuvent être intégrés dans une procédure plus générale et nous l'illustrons par la recherche de la frontière de Pareto dans un problème de clustering bi-critères sous contraintes utilisateur.Dans la tache de clustering, le critère de diamètre vise à trouver des clusters compacts mais soure souvent de l'eet de découpage[Cormack 1971], i.e. des objets assez similaires peuvent être classés dans des classes diérentes pour réduire le diamètre. Par contre, le critère de la marge peut sourir de l'eet de chaîne[Johnson 1967], i.e. des objets très diérents peuvent se trouver dans une même classe. Pour éviter ces problèmes, on peut considérer ces critères en même temps. Nous présentons dans ce chapitre la tâche de clustering bi-critères, à savoir maximiser la marge entre clusters et minimiser le diamètre maximal des clusters. Pour résoudre ce problème, une approche générale est de chercher la frontière de Pareto. Nos modèles orent à l'utilisateur le choix du critère d'optimisation et la capacité d'intégrer diérents types de contraintes utilisateur. Nous montrons que la généricité et la déclarativité de nos modèles permettent de l'utiliser directement pour traiter le clustering bi-critères de marge-diamètre et sous contraintes utilisateur. Ceci est réalisé par des optimisations mono-critère successives utilisant des contraintes utilisateur. Des expérimentations sur des bases de données classiques montrent que l'utilisation de notre modèle est plus ecace que les algorithmes exacts existants.

	Chapter 7
	Bi-criteria constrained clustering by
	CP
	Dans ce chapitre

Table 7 .

 7 1: Comparison of performance with bi-criteria as well as an extension with Large Neighborhood Search is proposed in [?]. This constraint

	7.4. Experimentation

  Classication non supervisée sous contrainte utilisateurs par la programmation par contraintes Résumé : La classication non supervisée, souvent appelée par le terme anglais de clustering, est une tâche importante en Fouille de Données. Depuis une dizaine d'années, la classication non supervisée a été étendue pour intégrer des contraintes utilisateur permettant de modéliser des connaissances préalables dans le processus de clustering. Diérents types de contraintes utilisateur peuvent être considérés, des contraintes pouvant porter soit sur les clusters, soit sur les instances. Dans cette thèse, nous étudions le cadre de la Programmation par Contraintes (PPC) pour modéliser les tâches de clustering sous contraintes utilisateur. Utiliser la PPC a deux avantages principaux : la déclarativité, qui permet d'intégrer aisément des contraintes utilisateur et la capacité de trouver une solution optimale qui satisfait toutes les contraintes (s'il en existe). Nous proposons deux modèles basés sur la PPC pour le clustering sous contraintes utilisateur. Les modèles sont généraux et exibles, ils permettent d'intégrer des contraintes d'instances must-link et cannot-link et diérents types de contraintes sur les clusters. Ils orent également à l'utilisateur le choix entre diérents critères d'optimisation. An d'améliorer l'ecacité, divers aspects sont étudiés. Les expérimentations sur des bases de données classiques et variées montrent qu'ils sont compétitifs par rapport aux approches exactes existantes. Nous montrons que nos modèles peuvent être intégrés dans une procédure plus générale et nous l'illustrons par la recherche de la frontière de Pareto dans un problème de clustering bi-critère sous contraintes utilisateur.

For a discrete variable,[i, j] denotes the set of integers from i to j.

[1, n] denotes the set of integers from 1 to n.

until no change of cluster centroid;To overcome the disadvantage of COP-Kmeans algorithm, an approach is to try not to satisfy all the constraints, but most of them. The CVQE (constrained vector quantization error) algorithm[Davidson & Ravi 2005b] extends the K-means algorithm to incorporate

Chapter 5. An initial CP Model for Constrained Clustering Figure 5.5: Example of search strategy for the criterion of diameter and split

Remerciements

Chapter 5. An initial CP Model for Constrained Clustering

The optimal diameter is guaranteed with an optimization constraint: minimize D.

(5.15)

The constraint (5.15) is activated each time a solution ∆ is found. The diameter maximum of the solution, denoted by D ∆ , is obtained by the value D.lb, then a new constraint is added to ensure that the next solution (if it exists) must have a better maximum diameter:

(5.16)

Modeling the split criterion

The criterion of split is encoded with the oat variable S. At the beginning, the domain of S is not a set of values, but an interval:

where d min and d max are the minimum and the maximum dissimilarity between any two points.

For the split criterion, a similar requirement is put on each pair of points: if the two points are in dierent clusters, then the dissimilarity between them must be equal or superior to the split:

It is equivalent with the constraint: if the dissimilarity between any two points is inferior to the split, then they must be in the same cluster.

Similar to the criterion of diameter, for each pair of points i, j, two auxiliary boolean variables b ij1 and b ij2 are generated and three constraints are posted as follows:

The propagation of these constraints is as follows: (5.21). The reied constraint (5.20) then propagates the binary constraint

), there is no propagation.

• If b ij2 = 1 or (G i = G j ) holds, there is no propagation.

• If b ij2 = 0 or (G i = G j ) holds, then the variable b ij1 is instantiated to 0, the constraint d ij ≥ S is propagated. The lower bound of the domain of S is modied:

These propagations ensure that the lower bound of S is the minimum split.

Criterion Modeling 63

The optimal split is guaranteed with an optimization constraint: minimize S.

With this optimization constraint, each time a solution ∆ is found, the minimum split S ∆ is obtained by the value of S.lb, then a new constraint is added:

S > S ∆ .

Modeling the criterion of Within Cluster Sum of Dissimilarities

The Within Cluster Sum of Dissimilarities (WCSD) is measured as the sum of dissimilarities between any two points in the same cluster, where the dissimilarities is usually dened as the squared Euclidean distance. To represent the relation between points and WCSD, we need a constraint wcsd(G,V,d), which ensures:

( 5.22) where

Using predened constraints, the constraint (5.22) for the WCSD criterion can be implemented by the following constraints:

However, these constraints, while considered independently, do not oer enough propagation. For example, with k = 2, given 4 points from 1 to 4 and a partial assignment where G 1 = 1 and G 2 = G 3 = 2 as in Figure 5.3 (the number on each edge {i, j} represents the value d ij ). We have three instantiated boolean variables:

The remaining boolean variables are uninstantiated:

Let us assume that a solution with V = 5 was found. With the branch-and-bound search, this solution sets the upper bound of variable V to 5. A new constraint is added:

As b 12 = b 13 = 0, b 23 = 1, this constraint becomes:

Chapter 5. An initial CP Model for Constrained Clustering Algorithm 11: Filtering algorithm

So if this last value is greater than the actual upper bound of V , v is inconsistent. The complexity of this algorithm is O(n 2 + nk), since the domain of each G i is of size at most k. Since k ≤ n, the complexity is then O(n 2 ).

Search Strategy

A search strategy must be given to indicate to the solver the choice of variables and values. Variables will be chosen in the order I then G. This order represents the fact that the cluster representatives must be identied rst, then the solver tries to determine the assignment of points to clusters. Variables in I are chosen from I 1 to I k , that means identifying the representative of the rst cluster till the last cluster. For the choice of values for each I c , since the representative is the smallest index in each cluster, the values are enumerated in the ascending order. The next step is branching on variables in G. The branching on uninstantiated variables in G nds a variable G i and a value c in the domain of G[i] and makes two alternatives:

For the criterion of maximum diameter and minimum split Variables in G are chosen in the ascending order of the size of their remaining domain. For the choice of values for each G i , the closest representative is enumerated rst.

Example 24 Given a dataset of 4 points and a clustering task dividing points into 2 clusters. Figure 5.5 expresses the search tree by using the explained strategy. Node 0 presents the initial domain of variables at the beginning of the search. As I 1 is instantiated, the variable I 2 is branched by assigning I 2 = 2 rst. When I 2 = 2, the variable G 2 is instantiated to 2 and the domain of others variables in G is restricted to {1, 2}. The variable G 3 is chosen for the next branching. As the representative point I 2 is actually closer to the point 3, G 3 is assigned rst to 2, as expressed in Node 2. In Node 19, as soon as the second representative point is instantiated I 2 = 4, the variables G 2 , G 3 and G 4 are all instantiated, by the propagation of the partition constraints: the representative point is the point with the smallest index in the cluster.

For the criterion of WCSD A mixed strategy is used. Because an upper bound is necessary for the WCSD constraint to be eective, a greedy strategy is used rst to quickly nd a solution. In this step, G i and c are selected to make sure that the value of V increases as little as possible. The rst solution found in general has a good value of WCSD (which is close to the optimal value). After nding a rst solution, the search strategy is changed, following the fail-rst principle, which tends to cause the failure early. In this strategy, the branching tries to make alternatives on frontier points. The value c is always the representative of the closest group to point i. The variables in G with the smallest domain size are all considered and the variable that makes the most changes on V is chosen for branching.

Modeling the WCSD criterion

For the WCSD criterion, we use the same constraint wcsd(G, W, d), presented in section 5.4.3 of Chapter 5. This constraint relies only on variables in G, and more specically, only on the condition G i = G j . The ltering algorithm for this constraint is therefore mostly unchanged in this model. The only dierence is, in this model, the number of clusters k is not xed but k min ≤ k ≤ k max . Therefore, for calculating the low bound of WCSD, the number of clusters k is set as the maximum number of clusters: k = max(∪ i∈[1,n] Dom(G i )).

Search strategy

To improve the eciency, points are beforehand reordered following the strategy presented in section 5.5 in the previous Chapter. For the search strategy, branching is realized on the variables of G.

For the diameter and split criteria At each branching point, a variable G i with the smallest remaining domain is chosen. With this variable G i , all the values c in the domain of G i are examined and the number of the closest cluster to i is chosen. The dissimilarity between a point i and a cluster number c is dened as the maximal dissimilarity d ij where G j is assigned and G j = c. If the cluster number c is empty (there is no point j such that G j = c), the dissimilarity between i and the cluster c is null. This means that the