An optimization-based decomposition heuristic for the microaggregation problem

An Optimization-Based Decomposition Heuristic for the Microaggregation Problem Jordi Castro1(B), Claudio Gentile2, and Enric Spagnolo-Arrizabalaga1 1Department of Statistics and Operations Research, Universitat Polit`ecnica de Catalunya, Jordi Girona 1–3, 08034 Barcelona, Catalonia [email protected],[email protected] 2Istituto di Analisi dei Sistemi ed Informatica “A. Ruberti”, Consiglio Nazionale delle Ricerche, Rome, Italy Abstract. Given a set of points, the microaggregation problem aims to find a clustering with a minimum sum of squared errors (SSE), where the cardinality of each cluster is greater than or equal to k.Pointsinthe cluster are replaced by the cluster centroid, thus satisfying k-anonymity. Microaggregation is considered one of the most effective techniques for numerical microdata protection. Traditionally, non-optimal solutions to the microaggregation problem are obtained by heuristic approaches. Recently, the authors of this paper presented a mixed integer linear optimization (MILO) approach based on column generation for computing tight solutions and lower bounds to the microaggregation problem. However, MILO can be computationally expensive for large datasets. In this work we present a new heuristic that combines three blocks: (1) a decomposition of the dataset into subsets, (2) the MILO column generation algorithm applied to each dataset in order to obtain a valid microaggregation, and (3) a local search improvement algorithm to get the final clustering. Preliminary computational results show that this approach was able to provide (and even improve upon) some of the best solutions (i.e., of smallest SSE) reported in the literature for the Tarragona and Census datasets, and k∈{3,5,10}. Keywords: Statistical disclosure control ·Microdata · Microaggregation problem ·Mixed integer linear optimization · Column generation ·Local search ·Heuristics 1 Introduction A microdata file of pindividuals (people, companies, etc.) and dvariables (or attributes) is, in practice, a matrix A∈Rp×dwhose element aij provides the value of attribute jfor individual i, and whose row aigives the dattributes for Supported by grant MCIU/AEI/FEDER RTI2018-097580-B-I00. c Springer Nature Switzerland AG 2022 J. Domingo-Ferrer and M. Laurent (Eds.): PSD 2022, LNCS 13463, pp. 3–14, 2022. https://doi.org/10.1007/978-3-031-13945-1_1 4J.Castroetal. individual i. Formally, a microdata file is a mapping M:S⊆P→V1×...×Vt, where Pis a population, Sis a sample of the population and Viis the domain of the attribute i∈{1,...,d}. Microdata files must be protected before being released; otherwise, confidential individual information would be jeopardized. Microaggregation [5,6]is a statistical disclosure control technique, mainly for numeric variables, which is related with k-anonymity [20]. The goal of microaggregation is to modify the values of the variables such that the released microdata satisfies k-anonymity. Therefore, it first partitions the individuals (or points in Rd) into subsets of size at least k, called clusters, and it then replaces each point in the cluster with the centroid of the cluster in order to minimize the loss of information, called spread. In practical cases, the value of kis relatively small (e.g., 3 ≤k≤10, see [6]). A widely used measure for evaluating the spread is the sum of squared errors (SSE)[6]: SSE = q  i=1 ni  j=1 (aij−ai)T(aij−ai),(1) where qdenotes the number of clusters, nithe size of cluster Ci={aij,j = 1,...,n i},andai=1 nini j=1 aijits centroid, for i=1,...,q. An equivalent measure that is also widely used in the literature is the information loss (IL), which is defined as IL =SSE SST ·100,(2) where SST is the total sum of squared errors for all the points, that is: SST = p  i=1 (ai−¯a)(ai−¯a) where ¯a=p i=1 ai p.(3) IL always takes values within the range [0,100]; the smaller the IL, the better the clustering. From now on, we will denote as feasible clustering a partition into clusters of size at least k. Finding the partition that minimizes IL (or SSE) and satisfies the cardinality requirement ni≥kfor i=1,...,q is a difficult combinatorial optimization problem when d>1 (multivariate data), which is known to be NP-hard [15]. For univariate data (that is, d= 1)—which in practice are the exception— microaggregation can be solved in polynomial time using the algorithm of [11], which is based on the shortest path problem. Microaggregation differs from other related clustering problems, such as kmedians or k-means [10], specifically in that it imposes a lower bound kon the cardinality of each cluster, but no fixed number of clusters. On the other hand, kin k-medians and k-means fixes the number of clusters, while imposing no constraint on the cardinality of each cluster. There exist various papers on heuristic algorithms for feasible solutions to multivariate microaggregation with reasonable IL. Heuristics like minimum distance to average (MDAV) [6,8]andvariable minimum distance to average Decomposition Heuristic for the Microaggregation Problem 5 (VMDAV) [19] sequentially build groups of fixed (MDAV) or variable (VMDAV) size based on considering the distances of the points to their centroid. Other approaches first order the multivariate points and apply the polynomial time algorithm of [11] to the ordered set of points, such as in [7], which used several fast ordering algorithms based on paths in the graph that is associated with the set of points, whereas [14] used (slower) Hamiltonian paths (which involve the solution of a traveling salesman problem). The heuristic of [16] also sequentially builds a set of clusters attempting to locally minimize IL. Other approaches, such as those of [2,13], are based on refining the solutions previously provided by another heuristic. To our knowledge, the only two papers in the literature to apply optimization techniques to microaggregation and formulate it as a combinatorial optimization problem are [1] and [4]. Both of them apply a column generation algorithm inspired by the work in [9]. Those optimization approaches solve the linear relaxation of the integer microaggregation problem, thus computing a (usually tight) lower bound for the problem. They also provide a (usually very good) upper bound solution with an IL that is smaller than the ones reported by other heuristics. Note that having a lower bound of the optimal solution is instrumental in order to know how good are the (upper bound) solutions computed by heuristics, even to perform fair comparisons between them. For instance, the heuristic introduced in [17] reported IL values below the certified lower bound— thus, not possible—, which clearly indicates that the values of the datasets used in that paper were different than those in the rest of the literature (likely due to some sort of normalization of attributes). The downside of those optimization based techniques is that, when the dataset is large, the column generation may involve a large number of iterations, thus making it computationally very expensive. The main difference between the approaches in [1] and [4] is that the pricing problem of the former involves a nonlinear integer problem while the latter requires a simpler linear integer problem. In practice this means that the pricing subproblem in [1] can be tackled only by means of complete enumeration and only for small values of k, while [4] theoretically offers more flexibility and can deal with larger values of k. Since optimization-based methods can be inefficient for large datasets but can provide high quality solutions in reasonable time for microdata with a small number of points, this work suggests a new approach consisting of first partitioning the set of points, and then applying an optimization approach to each smaller subset. The initial partitioning of the dataset is done according to a feasible clustering previously computed with the MDAV/VMDAV heuristics. Additionally, a local search improvement heuristic is also applied twice: first, to the solution provided by the MDAV/VMDAV heuristics prior to the partitioning; and second, to the final solution provided by the optimization approach. This short paper is organized as follows. Section 2outlines the optimizationbased decomposition heuristic. Sections 3and 4outline the two main building blocks of the heuristic: the local search improvement algorithm and the mixed integer linear optimization method based on column generation in [4]. Section 5 6J.Castroetal. shows the preliminary results from this approach with the standard Tarragona and Census datasets used in the literature. 2 The Decomposition Heuristic The decomposition heuristic comprises the following steps: Input: Microdata matrix A0∈Rp×d, minimum cluster cardinality k, upper bound of the number of subsets sin which the microdata will be decomposed. 1. Standardize attributes/columns of A0, obtaining matrix A∈Rp×d. Compute the squared Euclidean distance matrix D∈Rp×p, where Dij = (ai−aj)(ai−aj), which is to be used in the remaining steps. 2. Apply MDAV and VMDAV microaggregation heuristics using D.LetC= {C1,...,Cq}be the best of the two feasible clusterings provided by MDAV and VMDAV (that is, the one with smallest IL). Here, qrepresents the number of clusters, and Cithe set of points in cluster i. 3. Apply the local search improvement algorithm (described in Sect. 3)to C, obtaining the updated clustering C={C1,...,Cq}. The updated clustering Chas the same number of clusters qas C, but the points in subsets Ciand Cican be different. 4. Partition the microdata and distance matrices Aand Din s≤ssubsets Si,i=1,...,s , according to the clustering C. For this purpose we compute ¯p=round(p/s), the minimum number of points in each subset of the partition, and build each subset Siby sequentially adding points of clusters Cj,j =1,...,q until the cardinality ¯pis reached. 5. Apply the mixed integer linear optimization method based on column generation in [4] to each microaggregation subproblem defined by ASi and DSi,i=1,...,s . Obtain feasible clustering Oifor points in Si, i=1,...,s . 6. Perform the union of clusterings O=O1∪ ··· ∪ O s.Ois a feasible clustering for the original microdata A. 7. Finally, once again apply the local search improvement algorithm from Sect. 3to Oin order to obtain the final microaggregation O. Return: Clustering O. Step 3 of the algorithm can be skipped, thus obtaining the partition in Step 4 with the clustering Cfrom Step 2. However, we have observed that better results are generally obtained if the local search improvement heuristic is applied in both Steps 3 and 7, not only in Step 7. Indeed, if efficiency is a concern, it is possible to stop the whole procedure after Step 3, which thus returns cluster Cas a solution and, in general, significantly outperforms the solution obtained in Step 2. In this way, it is possible to avoid Step 5, which is usually computationally expensive. Note also that the clustering Cobtained in Step 3 is used only in Step 4 to decompose the microdata into subsets, but not as a starting solution for Decomposition Heuristic for the Microaggregation Problem 7 Fig. 1. Two-swapping local search improvement heuristic the optimization procedure in Step 5. Therefore, the clustering computed in Steps 5–6 by the optimization procedure might have a larger SSE than C.On the other hand, by not starting the optimization procedure in Step 5 with the solution Cwe have some chances to obtain a different and possibly better local solution. In the current implementation, Cis not used as a starting solution for the optimization algorithm. The larger the value of s, the faster the algorithm will be, since the minimum number of points ¯p=round(p/s) in each subset Si,i=1,...,s will be smaller, and therefore the optimization algorithm of [4] will be more efficient. However, the final IL (SSE) of the final clustering Oalso increases with s. Therefore parameter scan be used as a trade-off between efficiency and solution quality. In the next two sections, we outline the local search improvement heuristic used in Steps 3 and 7, as well as the mixed integer linear optimization method in Step 5. 3 The Local Search Improvement Heuristic Given a feasible clustering for the microaggregation problem, a local search algorithm tries to improve it by finding alternative solutions in a local neighborhood of the current solution. The local search considered in this work is a two-swapping procedure; in addition to its simplicity, it has proven to be very effective in practice. Briefly, the two-swapping heuristic performs a series of iterations, and at 8J.Castroetal. each iteration it finds the pair of points (i, j) located in different clusters that would most reduce the overall SSE if they were swapped. This operation is repeated until no improvement in SSE is detected. The cost per iteration of the heuristic is O(p2/2). Similar approaches have been used in other clustering techniques, such as in the partitioning around medoids algorithm for k-medoids [12]. The two-swapping algorithm implemented is shown in Fig. 1. 4 The Mixed Integer Linear Optimization Algorithm Based on Column Generation In this section we quickly outline the optimization method presented in [4]. Additional details can be found in that reference. The formulation of microaggregation as an optimization problem in [4]is based on the following property of the SSEhof cluster Ch={ahi,i=1,...,n h} (see [4, Prop. 3] for a proof): SSEh= nh  i=1 (ahi−ah)(ahi−ah) =1 2nh nh  i=1 nh  j=1 (ahi−ahj)(ahi−ahj)= 1 2nh nh  i=1 nh  j=1 Dhihj. (4) That is, for computing SSEh, we do not need the centroid of the cluster, but only the distances between the points in the cluster. From (4), defining binary variables xij ,i, j =1,...,p(which are 1 if points i and jare in the same cluster, 0 otherwise), then the microaggregation problem can be formulated as: min SSE 1 2 p  i=1 p j=1,j=iDij xij p j=1,j=ixij +1 (5a) s. to xir +xjr −xij ≤1i, j, r =1,...,p, i=j, r =j, i =r(5b) p  j=1,j=i xij ≥k−1i=1,...,p (5c) xij =xji,x ij ∈{0,1},i,j =1,...,p. (5d) Constraints (5b) are triangular inequalities, that is, if points iand r,andr and jare in the same cluster, then points iand jare also in the same cluster. Constraints (5c) guarantee that the cardinality of the cluster is at least k.The denominator in the objective function (5a) is the cardinality of the cluster that contains point i. Unfortunately, (5) is a difficult nonlinear and nonconvex integer optimization problem (see [4] for details). A more practical alternative is to consider a formulation inspired by the clique partitioning problem with minimum clique size of [9]. Defining as C∗={C ⊆ Decomposition Heuristic for the Microaggregation Problem 9 {1,...,p}:k≤|C|≤2k−1}the set of feasible clusters, the microaggregation problem can be formulated as: min  C∈C∗ wCxC s. to  C∈C∗:i∈C xC=1i∈{1,...,p} xC∈{0,1}C∈C ∗, (6) where xC= 1 means that feasible cluster Cappears in the microaggregation solution provided, and the constraints guarantee that all the points are covered by some feasible cluster. From (4), the cost wCof cluster Cin the objective function of (6)is wC=1 2|C| i∈C  j∈C Dij.(7) The number of feasible clusters in C∗—that is, the number of variables in the optimization problem (6)—is 2k−1 j=kp j, which can be huge. For instance, for p= 1000 and k=3wehave|C∗|=8.29 ·1012. However, the linear relaxation of (6) can be solved using a column generation technique, where the master problem is defined as (6) but it considers only a subset ¯ C⊆C ∗of the variables/clusters. The set ¯ Cis updated at each iteration of the column generation algorithm with new clusters, which are computed by a pricing subproblem. The pricing subproblem either detects that the current set ¯ Ccontains the optimal set of columns/clusters or, otherwise, it generates new candidate clusters with negative reduced costs. For small datasets and values of k, the pricing subproblem can be solved by complete enumeration; otherwise, an integer optimization model must be solved. The master problem requires the solution of a linear optimization problem. Both the linear and integer optimization problems were solved with CPLEX in this work. The solution of the linear relaxation of (6) provides a lower bound to the microaggregation problem (usually a tight lower bound). In addition, solving the master problem as an integer problem allows us to obtain a feasible solution to the microaggregation problem (usually of high quality). A thorough description of this procedure, and of the properties of the pricing subproblem, can be found in [4] and [18]. 5 Computational Results The algorithm in Sect. 2and the local search heuristic in Sect. 3have been implemented in C++. We used the code of [4] (also implemented in C++) for the solution of the mixed integer linear optimization approach based on column generation, as described in Sect. 4. A time limit of 3600 s was set for the solution of each subproblem with the column generation algorithm in Step 5 of the heuristic in Sect. 2. We tested the decomposition algorithm with the datasets Tarragona and Census, which are standard in the literature [3]. The results for Tarragona 10 J. Castro et al. and Census are shown, respectively, in Tables 1–2and 3–4. These tables show, for each instance and value k∈{3,5,10},theIL and CPU time for the main steps of the decomposition heuristic in Sect. 2. We also tried the different values s∈{40,20,10,5,2}for partitioning the dataset in Step 5. For Step 2, the tables also show which of the MDAV or VMDAV algorithms reported the best solution. The difference between Tables 1and 2(also between Tables 3and 4) is that the former show results with the Step 3, while in the latter this step was skipped. We remind the reader that the decomposition algorithm could be stopped after Step 3 with a feasible and generally good solution. However, the best IL values for each k, which are marked in boldface in the tables, are obtained after Step 7, although this means going through the usually more expensive Step 5. Table 1. Results for the Tarragona dataset considering Step 3 of the algorithm. The best IL for each kis marked in boldface. Instance kStep 2 Step 3 Step 5 Step 7 Alg IL CPU IL CPU sIL CPU IL CPU Tarragona 3VMDAV 15.85 0.6 15.00 1.9 40 14.96 0.2 14.85 0.1 20 14.83 0.4 14.81 0.1 10 14.68 0.9 14.65 0.3 514.57 2.8 14.53 0.4 214.51 4.2 14.50 0.2 Tarragona 5MDAV 22.46 0.5 20.74 5.2 40 20.74 1.0 20.73 0.5 20 20.74 7.9 20.73 0.3 10 20.47 103.8 20.40 1.5 520.32 1119.6 20.25 1.6 221.06 7207.8 20.46 4.0 Tarragona 10 MDAV 33.19 0.3 30.77 25.6 40 30.77 12119.7 30.77 0.1 20 33.03 61600.3 30.87 9.0 10 31.80 88899.1 30.56 13.4 533.32 9.8 30.79 17.4 233.20 22.8 30.80 17.2 From Tables 1,2,3and 4we conclude that: – In general, the smaller the kand larger the s, the faster the decomposition heuristic. In a few cases, however, smaller values of salso meant smaller CPU times; for instance, this is observed for Census, k= 10, and values s=20and s= 10. The explanation is that the maximum time limit was reached in some pricing subproblems in those runs, and therefore the CPU time increased with s. – In general, smaller ILs are obtained for smaller svalues, as expected. Decomposition Heuristic for the Microaggregation Problem 11 Table 2. Results for the Tarragona dataset without Step 3 of the algorithm. The best IL for each kis marked in boldface. Instance kStep 2 Step 5 Step 7 Alg IL CPU sIL CPU IL CPU Tarragona 3VMDAV 15.85 0.6 40 15.15 0.214.95 1.7 20 14.86 0.314.73 1.3 10 14.75 1.014.66 1.2 514.58 1.914.54 0.7 214.51 4.614.50 0.4 Tarragona 5MDAV 22.46 0.5 40 21.81 0.921.18 4.7 20 21.14 7.920.59 4.7 10 20.62 86.220.36 3.8 520.37 894.520.29 2.1 221.01 7208.120.77 6.0 Tarragona 10 MDAV 33.19 0.3 40 32.05 12597.230.82 19.9 20 33.18 60243.530.81 20.8 10 32.23 230644.630.55 20.9 533.20 13.030.84 20.9 233.21 21.230.83 20.6 Table 3. Results for the Census dataset considering Step 3 of the algorithm. The best IL for each kis marked in boldface. Instance kStep 2 Step 3 Step 5 Step 7 Alg IL CPU IL CPU sIL CPU IL CPU Census 3VMDAV 5.66 1.2 5.25 3.3 40 5.21 0.15.20 0.2 20 5.20 0.25.19 0.1 10 5.21 0.65.18 0.3 5 5.12 3.25.07 0.8 2 4.85 5.24.79 0.4 Census 5VMDAV 8.98 1.1 8.12 8.7 40 8.12 2.48.12 0.2 20 8.12 22.58.11 0.2 10 8.14 247.08.09 0.8 5 7.96 5334.67.84 2.0 2 9.36 7209.38.19 8.5 Census 10 VMDAV 14.04 0.8 12.36 38.4 40 12.36 8452.712.36 0.4 20 12.96 77221.112.32 7.0 10 12.84 37037.712.40 14.5 5 13.49 7310.412.63 28.8 2 14.45 26.312.46 37.1