Binary data arise in several situations in research since they can encode situations where the presence (1) or absence (0) of a characteristic is studied: species present/absent, alive/dead in health sciences, yes/no in social or decision sciences, and so on. For instance, in Ecology [14], it is usual to divide an area in sectors and record the presence or absence of certain vegetable species. In the study of gene expression data, several molecular techniques encode data as binary matrices [4]. In pattern recognition, images may also be coded as 0/1 data indicating the presence or absence of a feature. In Sociology [2], Health Sciences [18], Economics [10] . . . binary data are analyzed.

Several methods have been used for clustering binary data. For instance, there are partitioning methods such as dynamical clusters, which is an adaptation of Forgy’s k-means [6], based on a representation of clusters by a kernel and iterations of two steps: allocating objects to the nearest kernel and recalculating the kernels. A variant of this k-means method is PAM [11], partitioning around medoids, where kernels are 0/1 median vectors and L1 dissimilarity is used. Another variant is for kernels selected as the object in the class that minimizes the sum of dissimilarities to the rest of objects in the class; this last version is what we will call k-means for 0/1 data in this article. These methods find local minima of the criterion to be minimized, since they are based on local search procedures.

There are also hierarchical methods used for clustering binary data [6], using an appropriate dissimilarity for binary data (such as Jaccard, for instance). These methods find also local minima since they are based on a greedy strategy.

To overcome the local optima problem, optimization strategies have been used. We have employed combinatorial optimization metaheuristics for clustering numerical data [16], [17]. In the present article, we will used some of these heuristics for binary data, whenever it is possible. When dealing with binary data it is necessary to adapt the criterion, since Huygens theorem and other theoretical results only hold in an Euclidean context.

The article is organized as follows. Section 2 presents the clustering problem and particularly some criteria and properties for the binary case. Section 3 contains the five combinatorial optimization metaheuristics employed here and the adaptation we made for clustering binary data. In Section 5 the results obtained are presented, and finally in Section 6 some concluding remarks are made.

Given a data set of binary vectors Ω = with and a number , we seek for a partition P = () of Ω such that elements in a class are more similar than elements of other classes.

We define a within inertia measure of the partition P as

where ) can be defined (among others) as

d being a binary dissimilarity index in Ω (for instance, Jaccard or ) the median vector of . In [15] and [13] there are studied some other criteria for clustering. These indexes satisfy the following monotonicity property.

Proposition 1 (see [13]) Let P = () and = () be partitions of Ω in K and K + 1 non-empty classes, and on 2satisfy the monotonicty property, since for all instances of the data, we have

Proofs can be found in [13] or by request upon the authors. As a consequence of proposition 1 the number of clusters must be predefined, since best clusterings with different number of classes cannot be compared. In [13] it is also proved that for and all optimal partitions have non empty classes.

Analogously to the continuous case, total inertia can be defined as:

and, thanks to the monotonicity, the between-classes inertia is defined as:

and it is always positive.

We have implemented clustering algorithms using well-known combinatorial optimization heuristics. Three of them are based on neighbors, simlated annealing (SA), threshold accepting (TA) and tabu search (TS), and two based on a population of solutions, genetic algorithm (GA) and ant colony optimization (AC).

state, a neighbor is defined by the single transfer of an object from its current class to a new one, chosen according to the corresponding heuristic rules. This will be the case in SA, TA and TS, and the mutation operation in GA.

Simulated annealing is a random search algorithm [12] that uses an external parameter of control called temperature, that controls the random acceptance of states worse than the previous one. It is employed the Metropolis rule for this acceptation: a new state generated from P is accepted if ∆W < 0, where ∆), otherwise, it can be accepted with probability exp(. It is well known [1] that under a Markov chain modeling simulated annealing has asymptotic convergence properties to the global optimum under some conditions, so it is expected that its use permits to avoid local minima. We found some simplification properties for ∆W, calculated in each iteration and useful for speeding up the algorithm.

SA parameters are: the initial acceptation rate, L length of Markov chains, factor reduction for such that and stop criterion.

TA was proposed by [5] and can be seen as a particular case of SA, with a different rule for acceptation. Movements that produce an improvement for the objective function in a neighborhood are accepted and movements that worsen it are accepted if they fall into a threshold that is positive and decreases in time. Clearly, the acceptation rule is deterministic, not stochastic.

TA parameters are the initial threshold, the factor reduction for threshold T h, the maximum number of iterations and the stop criterion.

TS [7] handles a tabu list T of length |T | with solutions or codes of solutions, that are forbidden to be attained for a certain number of iterations. In each step, current state moves to the best neighbor outside T . In our partitioning problem, T stores a code of the transfers that define the neighbors, in this case the indicator function of cluster that contained the object that changed of class during the transfer; that is, all objects that were together in the previous state, are forbidden to be together for |T | iterations.

TS parameters are |T |, maximum number of iterations and sampling size of neighborhoods.

GA handles a population of solutions, which are chromosomes representing partitions. A chromosome is an n-vector in an alphabet of K numbers indicating the presence/absence of the i-th object in the corresponding class. The fitness function is defined as f(P) = . In the genetic algorithm iterations, chromo- somes are kept using a random wheel roulette with elitism, with a probability proportional to f. For crossover between two parents selected at random (with a uniform distribution), the dominant parent is the one with the greatest fit-ness; a class in it is selected uniformly at random and this class is copied in the other parent, generating a new son. Mutation is the classical one: an object is selected randomly, and it is transferred to a new class.

Parameters for GA are the number M of partitions (population size), probabilities of crossover and mutation, maximum number of iterations. Iterations stop when the fitness variance of the population is less than .

In AC [3], each ant m is associated with a partition P, which is modified during the iterations. Given an object , the probability of transfering another object to the same class as in iteration t is defined as

M being the number of ants, weights, and an evaporation parameter.

AC parameters to calibrate are , size of population M, precision and the maximum number of iterations.

We performed a Monte Carlo-type experiment, generating 16 data tables controlling the following factors: n, the number of objects (with levels 120 and 1200); K the number of clusters (levels 3 and 5); , the cardinality of the clusters (equal cardinalities and one big cluster with 50% of the objects and the rest of objects distributed equally); and , the probability of the Bernoulli distribution, with levels of well separated clusters (= 0.1, 0.5, 0.9 for K = 3 and = 0.05, 0.25, 0.5, 0.75, 0.95 for K = 5) or fuzzy separated clusters (= 0.3, 0.5, 0.7 for K = 3 and = 0.2, 0.35, 0.5, 0.65, 0.8 for K = 5). Table 1 presents the characteristics of each of the 16 data tables generated in the experiment, including the factors and the levels.

Table 1: Characteristics of data tables generated, 4 factors with 2 leves each. n: number of objects, K: number of clusters, : cardinality of cluster probability of membership to

We have compared the results obtained with the five metaheuristics (SA, TA, TS, GA, AC) with two classical methods: PAM (partitioning around medoids) when using the dissimilarity index, and hierarchical clustering (HC) using average linkage. For each heuristic, a parameter calibration was performed, exploring different values for each parameter. After this calibration, parameters selected for this article were:

• Simulated annealing: = 0.95, L = 50; = 0.91, = 0.01.

• Threshold accepting: = 100, maxiter = 50, = 0= 0.01.

• Tabu search: maxiter = 150, |T | = 5, s = 0.1|N(P)|.

• Genetic algorithm: = 0= 0.8, M = 20, maxiter = 500, = 0.01.

• Ant colony: = 0= 0= 0.5, M = 10, maxiter = 500, = 0.01.

In Table 2 we report, for a multistart of size 100, the best value of W obtained so far by any method (noted ), the mean value of W for each method (noted W) and the attraction rate or percentage of times that was obtained (up to a relative error of 0.05).

Table 2: Results summary with for a multistart of size 100. is the best value obtained by any method, the attraction rate of for each method, and W the mean value for each method

Generally speaking, with simulated annealing we obtain good results, although PAM obtains good results in some cases. Threshold accepting sometimes reaches the optimum and tabu search only in two cases. Population based heuristics did not get good results with our implementation. It is worth noting that hierarchical clustering never obtained the optimum. Even if we do not report running times, SA is fast enough to be competitive. The main drawback of using heuristics, is tuning the parameters, though SA may have a standard choice.

This research was partially supported by the University of Costa Rica (CIMPA project 821-B1-122 and the Graduate Program in Mathematics) and the Costa Rica Institute of Technology. The supports are gratefully acknowledged.

[1] Aarts, E.; Korst, J.: Simulated Annealing and Boltzmann Machines. John Wiley & Sons, Chichester (1990)

[2] Borkulo, C.D. van, Borsboom, D., Epskamp, S., Blanken, T.F., Boschloo, L., Schoevers, R.A., Waldorp, L.J.: A new method for constructing networks from binary data. Scient. Rep. 4 (2014) doi: 10.1038/srep05918

[3] Bonabeau, E., Dorigo, M., Therauluz, G.: Swarm Intelligence. From Nat- ural to Artificial Systems. Oxford University Press, New York (1999)

[4] Demey, J.R., Vicente-Villardn, J.L., Galindo-Villardn, M.P., Zambrano, A.Y.: Identifying molecular markers associated with classification of genotypes by external logistic biplots. Bioinform. 24(24), 2832–2838 (2008)

[5] Dueck, G.; Scheuer, T.: Threshold accepting: A general purpose opti- mization algorithm appearing superior to simulated annealing, J. Comput. Phys. 90(1), 161–175 (1990)

[6] Everitt, B.S.: Cluster Analysis. Edward Arnold, London (1993)

[7] Glover, F.: Tabu search Part I. ORSA J. on Comput., 1, 190–206 (1989).

[8] Goldberg, D.E.: Genetic Algorithms in Search Optimization and Machine Learning. Addison-Wesley, Reading (1989)

[9] Jajuga, K.: A clustering method based on the -norm. Comput. Stat. & Data Analysis 5(4), 357–371 (1987) doi: 10.1016/0167-9473(87)90058-2

[10] Jeliazkov, I., Rahman, M.A.: Binary and ordinal data analysis in Eco- nomics: Modeling and estimation. In: Yang, X.S. (ed.) Mathematical Modeling with Multidisciplinary Applications, pp. 1–31. John Wiley & Sons, New York (2012)

[11] Kaufman, L., Roosseuw, P.: Finding Groups in Data. An Introduction to Cluster Analysis. John Wiley & Sons, New York (2005)

[12] Kirkpatrick, S., Gelatt, D., Vecchi, M.P.: Optimization by simulated an- nealing. Science, 220, 671–680 (1983)

[13] Piza, E., Trejos, J., Murillo, A.: Clustering with non-Euclidean distances using combinatorial optimisation techniques. In: Blasius, J., Hox, J., de Leeuw, E., Schmidt, P. (eds.) Social Science Methodology in the New Millennium,paper number P090504. Leske + Budrich, Darmstadt (2002)

[14] Salas-Eljatiba, C., Fuentes-Ramirez, A., Gregoire, T.G., Altamirano, A., Yaitula, V.: A study on the effects of unbalanced data when fitting logistic regression models in ecology. Ecological Indicators, 85, 502–508 (2018)

[15] Sp¨ath, H.: Cluster Dissection and Analysis. Theory, Fortran programs, Examples. Ellis Horwood, Chichester (1985)

[16] Trejos, J.; Murillo, A.; Piza, E.: Global stochastic optimization techniques applied to partitioning. In: Rizzi, A., Vichi, M., Bock, H.-H. (eds.) Advances in Data Science and Classification, pp. 185–190. Springer, Berlin (1998)

[17] Trejos, J., Villalobos, M., Murillo, A., Chavarr´ıa, J., Fallas, J.J.: Eval- uation of optimization metaheuristics in clustering. In: Travieso, C.M., Arroyo, J., Ram´ırez, M. (eds.) IEEE International Work-Conference on Bioinspired Intelligence, pp. 154–161. IEEE, Liberia (2014) doi: 10.1109/IWOBI.2014.6913956

[18] Zhang, H., Singer, B.: Recursive Partitioning in the Health Sciences. Springer, New York (1999)