Enhancing Scalability of Metric Differential Privacy via Secret Dataset Partitioning and Benders Decomposition

1 Introduction

Among the array of data privacy protection mechanisms, data perturbation has emerged as a widely employed technique to protect individuals' information. Data perturbation involves intentionally introducing noise into a database, rendering individual information unreadable by unauthorized users even in the event of data breaches on the server side. Particularly, Differential Privacy (DP) [Dwork et al., 2006] has become a paradigm of choice for data perturbation due to its strong and theoretically provable privacy guarantees.

In its original definition, DP requires data perturbation to maintain a uniform indistinguishability level for the queries of the neighboring databases [Dwork et al., 2006]. Specifically, the classification of "neighboring databases" relies on their Hamming distance, i.e., two databases are considered "neighbors" if they differ by at most one record. However, this definition limits DP's applicability in the domains where data similarity is evaluated by other distance metrics, and also varying degrees of indistinguishability are necessitated by nonbinary distance values among neighbors. Recognizing these constraints, DP has been extended to metric DP (mDP) [Imola et al., 2022], which generalizes the classification of "neighboring databases" by accounting for the diverse underlying distance metric spaces.

Related work of mDP. mDP originated in the domain of geo-location privacy protection [Andr´es et al., 2013], requiring "geo-indistinguishability" for each pair of locations with the Euclidean distance lower than a predetermined threshold. In other words, mDP defines "neighboring locations" based on Euclidean distance, diverging from the original DP which relies on Hamming distance. Over time, mDP has been explored across a spectrum of metric choices, including Manhattan distance [Chatzikokolakis et al., 2013], Hyperbolic distance [Feyisetan et al., 2019], Haversine distance [Pappachan et al., 2023], Word Mover's distance [Fernandes et al., 2019], and others [Feyisetan and Kasiviswanathan, 2021].

Compared to the conventional DP, the optimization of mDP poses additional challenges due to the diverse sensitivity of utility loss to data perturbation in general distance metric spaces [Qiu et al., 2022a]. A commonly employed approach is to discretize both the secret dataset and the corresponding perturbed data domain into finite sets [Imola et al., 2022], which allows for explicit measurement of the utility loss caused by each perturbation choice. In this case, the problem of optimizing the perturbation distribution of each secret record can be formulated as a linear programming (LP) problem [Fawaz and Shin, 2014], of which the objective is to minimize the expected utility loss caused by data perturbation while satisfying the mDP constraints for each pair of neighboring records.

However, the LP formulation of mDP typically requires O(N 2 ) decision variables and O(N 2K) linear constraints [Bordenabe et al., 2014], where N and K are the size of the secret dataset and perturbed dataset, respectively. The high complexity of LP makes it hard to apply in a large-scale mDP. As illustrated by Fig. 1, most current LP-based mDP works have to limit their applications to small-scale secret datasets CCS 2012 [Shokri et al., 2012], CCS 2014 [Fawaz and Shin, 2014], ICDM 2016 [Wang et al., 2016], WWW 2017 [Wang et al., 2017], NDSS 2017 [Yu et al., 2017], ICDCS 2019 [Qiu and Squicciarini, 2019], CIKM [Qiu et al., 2020], TMC 2022 [Qiu et al., 2022a], SIGSPATIAL 2022 [Qiu et al., 2022b], UAI 2022 [Imola et al., 2022], EDBT 2023 [Pappachan et al., 2023].

(i.e., up to 100 records in the secret dataset [Chatzikokolakis et al., 2015]). While some recent efforts [Qiu et al., 2020; Qiu et al., 2022a] have managed to extend the size of secret datasets to around 300 records using Dantzig-Wolfe decomposition and column generation algorithms, their focus has been on specific LP formulations where the mDP constraints are imposed for all pairs of secret records, simplifying the initialization of column generation. This specific LP formulation, however, is hard to apply in general mDP problems where the constraints are required only for neighboring records. Alternative approaches like [Imola et al., 2022] integrate LP with the exponential mechanism to enhance scalability, which, however, comes at the cost of sacrificing the optimality of mDP.

Our main contributions. Our primary aim is to enhance the scalability of the LP-based mDP by introducing a new computation framework. We first construct an mDP graph to describe the mutual mDP constraints among secret records and then partition the entire secret dataset into well-balanced subsets based on the mDP graph (Contribution 1). Building upon this partition, we then reformulate the LP and solve it using Benders Decomposition (Contribution 2), which consists of two stages: (1) a master program to manage the perturbation distribution calculations across subsets and (2) a set of subproblems, each deriving the perturbation distributions of records within a subset. The problems in both stages exhibit a relatively small scale, enabling efficient solutions. This efficiency facilitates the iterative derivation of a nearoptimal solution for the original LP through the communication between the two stages.

Finally, we test the performance of the new computation framework using geo-location data (in the road network and grid maps), text embeddings, and synthetic datasets, with a comparison of several benchmarks (Contribution 3). The experimental results demonstrate that our new computation framework can derive the optimal mDP for the secret dataset with 1,000 records, marking approximately 9 times improvement in scalability compared to the state-of-the-art methods for general mDP problems, as shown in Fig. 1.

2 Preliminaries

In this part, we introduce the preliminary knowledge of mDP in Section 2.1 and its LP computation framework in Section 2.2. Table 4 in [Qiu, 2024] lists the main notations used throughout this paper.

2.1 Metric DP

In general, a data perturbation mechanism can be represented as a probabilistic function Q: R → O, of which the domain R is users' secret dataset and the range O is the perturbed dataset. Given the underlying data features of R, a metric d : R2 ↦→ R is defined to measure the distance between any two records ri, rj ∈ R, where the distance between ri and rj is denoted by di,j .

Definition 2.1. (Neighboring records) Given the secret dataset R and its distance measure d, any pair of records ri, rj ∈ R are called neighboring records if their distance di,j ≤ η, where η > 0 is a pre-determined threshold.

In what follows, we use E = {(ri, rj) ∈ R2 |di,j ≤ η } to represent the whole set of neighboring records in R.

Definition 2.2. (Metric DP) For each pair of neighboring records, (ri, rj) ∈ E, ϵ-mDP ensures that the probability distributions of their perturbed data Q(ri) and Q(rj) are sufficiently close so that it is hard for an attacker to distinguish ri and rj according to the probability distributions of Q(ri) and Q(rj )

Here, ϵ > 0 is called the privacy budget, reflecting how much information can be disclosed from the perturbed data, i.e., lower ϵ implies a higher privacy level.

2.2 Linear Programing-based Approaches

Considering the computational challenges of optimizing a perturbation function Q defined on a continuous domain R and a continuous range O, it is common practice to discretize both R and O to finite sets [Imola et al., 2022]. In this case, the function Q can be represented as a stochastic perturbation matrix Z = {zi,k}(ri,ok)∈R×O , where each entry zi,k represents the probability of selecting ok ∈ O as the perturbed record given the real record ri ∈ R. As such, the mDP constraints formulated in Eq. (1) can be rewritten as the following linear constraints

Additionally, the sum probability of perturbed record ok ∈ O for each real record ri should be equal to 1 (the unit measure of probability theory [Stroock, 2010]), i.e.,

We let zi = [zi,1, ..., zi,K] (i = 1, ...,K) denote the record ri's perturbation vector, i.e., the probability distribution of its perturbed records. We let ci,k denote the utility loss caused by the perturbed record ok when the real record is ri. Then, the expected utility loss caused by the perturbation matrix Z can be represented by ∑Ni=1 ciz⊤i , where ci = [pici,1, ..., pici,K] and pi denotes the prior probability of the record being ri.

Consequently, the objective of the perturbation matrix optimization (PMO) problem is to minimize the expected data utility loss ∑Ni=1 ciz⊤i while satisfying the constraints of both mDP and the probability unit measure, which can be formulated as the following LP problem:

The decision variables in the LP problem in Eq. (4)-(7) are the perturbation matrix Z, including O(NK) decision variables (entries), constrained by O(N 2K) linear constraints. Such a high complexity makes this LP computation framework hard to apply in large-scale mDP applications [Pappachan et al., 2023].

3 Methodology

In this section, we present our approach to enhance the scalability of the mDP calculation via LP decomposition.

We first present the computational framework in Section 3.1, followed by the detailed descriptions of its two components: Secret dataset partitioning and Benders decomposition (BD). In our framework, the secret dataset partitioning precedes BD, but since the objective for partitioning the secret data is contingent on the complexity analysis of BD, we introduce BD in Section 3.2 and its time complexity analysis in Section 3.3 before delving into secret dataset partitioning in Section 3.4.

3.1 Computation Framework

In general, the decomposition of an optimization problem highly depends on how its decision variables are coupled by constraints [Palomar and Chiang, 2006]. A common strategy involves partitioning the original large-scale problem into a set of subproblems with smaller sizes. Within each subproblem, the decision variables are strongly linked by constraints, while the connections between decision variables across different subproblems are comparatively weaker.

In PMO, the decision variables are the perturbation vectors of secret records, zi = [zi,1, ..., zi,K] (i = 1, ...,K), linked by the mDP constraints (Eq. (2)). Here, we describe how the perturbation vectors of each pair of records are coupled by the mDP constraints in R by an undirected graph, called the mDP graph (Definition 3.1):

Definition 3.1. (mDP graph) The mDP graph G is defined as the ordered pair G = (R, E), where the node set R is the secret dataset and the edge set E is the whole set of neighboring records in R. The weight assigned to each edge (ri, rj) is the distance di,j between the two records.

As Fig. 2 shows, given the mDP graph G, our intuitive idea is to first partition G into M subgraphs: Gl = (R1, E1) , ..., GM = (RM, EM), such that the nodes within The mDP graph 𝒢 is partitioned to 𝒢ଵ, 𝒢ଶ,…, 𝒢𝑀

each subgraph are strongly connected by the mDP constraints (i.e., the perturbation vector of one node exerts a significant mDP constraint on the other), and the nodes across different subgraphs are weakly connected by the mDP constraints. Subsequently, a subproblem, denoted by Subl, can be formulated to determine the perturbation vectors of the secret records within each subgraph Gl (l = 1, ..., N). Here, we use zRl = {zi, |ri ∈ Rl} to denote the perturbation vectors of the set of records in Rl.

Next, we introduce Property 3.1 and Property 3.2 of the mDP graph to facilitate the computation of the PMO problem. The detailed proofs of the two properties can be found in Section B.1 and Section B.2 in [Qiu, 2024].

Property 3.1. (Independent computation of mDP components) Suppose that the mDP graph G has m components, denoted by Cl = (Rl, El) (l = 1, ...,m). Here, each component is a connected subgraph that is not part of any larger connected subgraph [Wilson, 1986]. Consider a subproblem Subl formulated to determine the optimal perturbation vectors for each Rl, where only the mDP constraints in El are considered. In this case, we can solve Sub1, ..., SubM independently, and their collective optimal solutions form the optimal solution for PMO.

Property 3.2. Given that two records ri and rj (not necessarily neighbors) are connected by a path in the mDP graph G, of which the shortest path distance (i.e., the sum weight of the edges in the path) is Di,j , then zi,k and zj,k are restricted by the following constraints:

In what follows, we restrict our attention to the case where the mDP graph G has a single component, considering the case of multiple components as straightforward to generalize according to Property 3.1. When G is a single-component graph, although dividing G facilitates the parallel computation of the major perturbation vectors within each subset, it remains imperative to jointly derive the perturbation vectors that are linked by mDP constraints across multiple subsets. For instance, in Fig. 2, the perturbation vector zi of ri in G1 must adhere to mDP constraints with the perturbation vector zj of rj in G2. Based on whether the records in each Rl are adjacent to records in other subsets by mDP constraints, we classify the records in each Rl into either "boundary records" or "internal records", as defined in Definition 3.2:

Definition 3.2. Each ri ∈ Rl is a boundary record if it has at least a neighbor rj in another subset Rn (n ̸= l) in the mDP graph G; otherwise, ri is an internal record.

We use Xl and Yl (Xl, Yl ⊆ Rl) to represent the internal record set and the boundary record set in Rl, respectively. We let zXl and zYl denote the perturbation vectors of Xl and Yl, i.e., zXl = {zi |ri ∈ Xl } and zYl = {zi |ri ∈ Yl }.

The PMO problem reformulation. After categorizing the records inRl into Xl and Yl (l = 1, ...,M), the original PMO formulated in Eq. (4)-(7) can be rewritten in a block ladder structure, as shown in Fig. 3:

1. The objective (OBJ) ∑Ni=1 cizi is rewritten as the sum of
2. ∑Ml=1 cYl zYl in the red block, representing the data utility loss of the boundary records in Y1, ..., YM, and
3. cXl zXl (l = 1, ...,M) in the blue blocks, representing the data utility loss of the internal records in Xl;
4. The constraints (CSTR) includes

constraints that connect the boundary perturbation vectors zY1 , ..., zYM across R1, ..., RM, and the unit measure constraints of each zYl (l = 1, ...,M),

• BYl zYl + AXl zXl ≥ bXl,Yl (l = 1, ...,M), including the mDP constraints between the perturbation vectors within Rl (including both Xl and Yl), and their unit measure constraints.

Such block ladder structure lends the reformulated PMO well to Benders decomposition (BD) [Rahmaniani et al., 2017]. Due to the limit of space, we list the detailed formulations of the coefficient matrices BY1,...,YM , BYl , AXl , and coefficient vectors bY1,...,YM , and bXl,Yl in Section A.2 in [Qiu, 2024].

3.2 Benders Decomposition

BD is composed of two stages, a master program (MP) to derive the perturbation vectors of all the boundary records, and a set of subproblems, labeled as Subl (l = 1, ...,M), to derive the perturbation vectors of the internal records in each subset Rl.

Stage 1: Master program. The MP derives the boundary records' perturbation vectors zY1 , ..., zYM and replaces the data utility loss of the internal records in each Xl by a single decision variable wl, i.e., wl = cXl zXl . The MP is formally formulated as the following LP problem

where each cut in H is a linear inequality of the decision variables zY1 , ..., zYM , w1, ...,wM. According to PMO's reformulation (in Fig. 3), wl is given by

Since the MP doesn't know the optimal values of zXl , instead of using Eq. (13), it "guesses" the value of wl based the cut set H. In the subsequent Stage 2, each Subl verifies whether the "guessed" value of wl is feasible and achieves the minimum data utility loss as defined in Eq. (13); if not, Subl proposes the addition of a new cut to be included in H, thereby guiding the MP to refine wl during the next iteration. Initial cut. According to the Property 3.2 of mDP constraints, the MP can initialize the cut set H by

where Di,j is the shortest path distance between ri and rj in the mDP graph G. In the following, we use {zY1 , ..., zYM ,w1, ...,wN} to represent the optimal solution of the MP.

Stage 2: Subproblems. After the MP derives its optimal solution {zY1 , ..., zYM ,w1, ...,wN} in Stage 1, each Subl validates whether wl has achieved the minimum data utility loss,

of which the dual problem can be formulated as the following LP problem:

There are three cases of the dual problem:

• C1: The optimal objective value is unbounded: By weak duality [Hillier, 2008], zYl does not satisfy AXl zXl ≥ bXl,Yl -BYl zYl for any zXl ≥ 0. Since the dual problem is unbounded, there exists an extreme ray ˜aXl s.t. A⊤Xl ˜aXl ≤ 0 and (bXl,Yl -BYl zYl )⊤ ˜aXl > 0. To ensure ˜aXl not to be an extreme ray in the next iteration, Subl suggests a new cut h (feasibility cut) to the MP:

• C2: The optimal objective value is bounded with the solution aXl : By weak duality, the optimal value of the dual

problem is equal to the optimal value of wl constrained on the choice of zYl . In this case, Subl checks whether wl < (bXl,Yl -BYl zYl )⊤ aXl +cYl zYl . If yes, then wl <minzXl ≥0 {cXl zXl |AXl zXl ≥ bXl,Yl -BYl zYl }, indicating that wl guessed by the MP is lower than the minimum data utility loss. Therefore, Subl suggests a new cut h : wl ≥ (bXl,Yl -BYl zYl )⊤ aXl +cYl zYl to the MP to improve wl in the next iteration.

C3: There is no feasible solution: By weak duality, the primal problem either has no feasible/unbounded solution. The algorithm terminates.

After adding the new cuts (from all the subproblems) to the cut set H, the BD moves to the next iteration by recalculating the MP and obtaining updated zYl and wl (l = 1, ...,M). As Stage 1 and Stage 2 are repeated over iterations, the MP collects more cuts from the subproblems, converging the solution zYl and wl to the optimal.

Proposition 3.3. (Upper and lower bounds of PMO's optimal) [Rahmaniani et al., 2017]

(1) The optimal solution of the MP (Eq. (9) - (12)) offers a lower bound of the optimal solution of the original PMO (Eq. (4)-(7)) (as the MP relaxes the constraints of PMO).

(2) The solution of the subproblems (Eq. (15)-(16)), if it exists, combined with the solution zYl of the MP, provides an upper bound of the PMO's optimal (since their solutions form a feasible solution of the original PMO).

Note that the optimal solution for PMOmust lie within the gap between the upper bound and the lower bound in Proposition 3.3. The smaller the gap between these two bounds, the closer the solution of the BD is approaching the optimal solution of PMO. Fig. 15-Fig. 18 (in Section D.5 in Appendix) gives some examples of how the upper and lower bounds change over iterations. Considering a prolonged convergence tail, we termine the algorithm when the gap between the best upper and lower bounds falls below a predetermined threshold ξ (e.g., ξ = 0.01 in our experiments).

3.3 Time Complexity Analysis

In the two-stage framework of BD, the number of decision variables in the MP in Stage 1 is K(∑Ml=1 |Yl|)+M, including K(∑Ml=1 |Yl|) variables (zY1, . . . , zYM ) to determine the perturbation vectors of the boundary records, and the M variables (w1, ...,wM) to "guess" the utility loss of subproblems. Each Subl in Stage 2 has K|Xl| decision variables to determine the perturbation vectors of the internal records zXl .

The MP and the subproblems are both formulated as LP problems, with time complexity depending on the number of decision variables [Cohen et al., 2019]. To facilitate analysis, we use a monotonically increasing function O(T(n)) to represent the time complexity of LP given the number of decision variables n. In each iteration, the MP in Stage 1 has O(T(K ∑Ml=1 |Yl|)) time complexity, and each subproblem l in Stage 2 (run in parallel) has O(T(K|Xl|)) time complexity. Since Stage 2 is terminated only after all subproblems are completed, its time complexity is O(maxl T(K|Xl|)) = O(T(Kmaxl |Xl|)) assuming all the subproblems are run in parallel. Hence, the total computation time of each iteration is given by O(T(K ∑Ml=1 |Yl|) + T(Kmaxl |Xl|)). Considering that both ∑Ml=1 |Yl| and maxl |Xl| is much smaller than the total number of records |R|, the time complexity of each iteration of BD is significantly lower than that of the original PMO. There are two questions remain:

1. How many iterations are needed to converge the solution to the optimal? In theory, if | ∪Ml=1 Yl| is finite, generalized BD (including LP formulations) ends within a finite number of iterations for any given ξ > 0 [Floudas, 2009]. To assess practical applicability, we will examine the convergence of our specific BD framework in Section 4.3 using different partitioning algorithms based on multiple datasets.

(2) How to optimize the dataset partitioning to minimize the computation time O(T(K ∑Ml=1 |Yl|)+T(Kmaxl |Xl|)) in each iteration? We introduce the detailed methods to address this problem in Section 3.4.

3.4 Secret Dataset Partitioning

According to the time complexity analysis of the BD framework in Section 3.3, we outline two primary objectives O-1 and O-2 for enhancing its computational efficiency when partitioning the secret dataset:

(O-1) Maintain strong mDP constraints within each subset and weak mDP constraints across different subsets. This not only facilitates relatively independent computation among subproblems but also reduces the number of decision variables (boundary records) in the MP.

(O-2) Balance the number of records in different subsets, to decrease the maximum time complexity of the subproblems.

Achieving the above two objectives is similar to solving the RatioCut problem [Hagen and Kahng, 1992], which seeks to partition a graph with minimal cuts while ensuring a wellbalanced size of subgraphs, a problem known to be NP-hard.

Considering the computational trajectability, we apply a Distance-Vector (DV)-based dataset partitioning algorithm. As the mDP constraints of two perturbation vectors depend on their records' distance, we embed each record ri by the distance vector di = [di,1, ..., di,N] to characterize its mDP relationship with all other records in the mDP graph. We consider that two records, ri and rj , exhibit a stronger connection when the Euclidean distance between their distance vectors,

∥di - dj∥2 = √∑Nl=1 (di,l - dj,l)2 , is lower, meaning that the two records are strongly coupled by the mDP constraints and also share similar mDP constraint relationship with other records. Accordingly, we partition the dataset to the subsets R1, . . . , RM using the distance vectors, which can be cast as the following k-means clustering formulation:

︸ ︷︷ ︸ Reflect the connection within Rl

where µl represents the centroid of the distance vectors in Rl, i.e., µl = ∑ri∈Rl di |Rl | (l = 1, . . . ,M). Note that the objective function in Eq. (18) attains a lower value when the records within each partition are strongly connected and the sizes of R1, ..., RM are well-balanced.

Besides the DV-based partitioning algorithm (labeled as "k-mean-DV" or "k-m-DV"), for comparison, we carried out three other benchmarks, record-based partitioning, adjacency matrix-based partitioning, and balanced spectral clustering [Hagen and Kahng, 1992], labeled as "k-m-rec", "k-madj", and "BSC", respectively. The details of the benchmarks are introduced in Section C in [Qiu, 2024].

4 Performance Evaluation

In this section, we test the performance of our approach using multiple datasets. We introduce the experiment settings in Section 4.1, and evaluate the performance of our method in terms of partitioning balance in Section 4.2 and computation efficiency in Section 4.3. More comprehensive experimental results (including the comparison with exponential mechanisms and the performance evaluation with different parameters, etc.) can be found in Section D in [Qiu, 2024] 1 .

4.1 Settings

Datasets. As mDP has been primarily used in text embeddings and geo-location data [Imola et al., 2022], we specifically choose text and geo-location datasets for the performance evaluation. Additionally, we assess our methods on a synthetic dataset.

(1) Geo-location dataset in the road network, composed of a set of "nodes" in the road network (retrieved by OpenStreetMap [OpenStreetMap, 2024]), including road intersections, forks, junctions where roads intersect with others, and points where the road changes direction. We selected the city Rome, Italy as the target region (specifically, the bounding area with coordinate (lat = 41.66, lon = 12.24) as the south-west corner, and coordinate (lat = 42.10, lon = 12.81) as the north-east corner). The target region encompasses a total of 43,160 discrete locations within Rome. In each experiment, we selected a relatively small target region in Rome, within which we sample 500 discrete locations in the road network. The distance metric between locations is defined by Haversine distance, i.e., the angular distance between two points on the surface of a sphere.

(2) Geo-location dataset in 10 grid maps, where each grid map is composed of 20×25 grid cells (the size of each cell is 1km×1km) in the target region, Rome, Italy. The center of each grid cell serves as a proxy for the cell's location. The distance metric between cells is defined as the Euclidean distance between their centers. In both datasets (1) and (2), we set η = 2.00km and ϵ = 10.0/km by default.

(3) Text dataset, which comprises 2,000 words. Each word is represented by a 300-dimensional vector using MATLAB's word2vec function [MATLAB, 2024b] (pre-trained for 1 million English words), capturing both the semantic meaning of the word and its contextual usage. The distance between words is measured by the Euclidean distance between their vectors. (4) Synthetic dataset, which comprises 2,000 3dimensional vectors, following the multivariate Gaussian distribution. The distance metric is defined as vectors' Euclidean

1 The MATLAB source code of our method is available at: https: //github.com/chenxiunt/MetricDP BendersDecomposition

Maximum size

Maximum size

distance. In both datasets (3)&(4), we set η = 2.0, ϵ = 10.0 by default.

Metrics: (1) Problem sizes of both subproblems and MPs. We measure the size of each decomposed problem by counting the number of perturbation vectors (each corresponds to a record) to derive. This metric reflects how effectively the dataset partitioning algorithms achieve balance in dividing the secret dataset. Note that the mDP graph of the MP may consist of multiple independent components. Hence, we measure the size of each component in the MP rather than the entire MP.

(2) Computation time to execute algorithms. The experiments are performed by a desktop with 13th Gen Intel Core i7 processor, 16 cores.

(3) Number of iterations for the BD to converge to the optiml solution. Considering the prolonged convergence tails of BD, we end the algorithm when the gap between its upper and lower bounds is lower than the threshold ξ = 0.01.

4.2 Secret Dataset Partitioning

Table 1 and Table 2 compare the sizes of subproblems and MP's components generated by the dataset partitioning algorithms, "k-m-DV", "k-m-rec", "k-m-adj", and "BSC". Each method undergoes testing on 10 instances for every dataset. The tables present the mean value±1.96×standard deviation of the average/maximum size of subproblems and MP components across different instances. The detailed visual representation of the dataset partitioning in different datasets is given in Fig. 11-Fig. 14 in [Qiu, 2024].

The tables show that, in all four datasets, k-m-DV partitions the PMO (containing 500 secret records) into sufficiently small subproblems and MP components, with the mean maximum size not exceeding 50. This decomposition facilitates efficient resolution of the individual problems using classical LP solvers such as the dual-simplex method [Hillier, 2008]. In contrast to k-m-DV, the size variation in the decomposed problems is higher for the other three methods. Specifically, the maximum problem size of the decomposed problems for k-m-adj and BSC is higher than 100, indicating a challenge in solving these larger problems efficiently. Both k-m-adj and BSC are less effective in balancing the size of decomposed problems, as they focus on the mDP constraint between neighboring records without considering overall constraint distribution. While k-m-rec demonstrates comparable performance to k-m-DV when applied to the geo-location data and synthetic data, it falls short in achieving balanced partitioning for text data, due to its limited ability to capture mDP constraints in high-dimensional (300-dimensional) spaces.

4.3 Computational Efficiency

Table 3 compares the computation time of the BD framework using the four dataset partitioning algorithms, "k-m-DV", "km-rec", "k-m-adj", and "BSC", and two classic LP solvers, dual-simplex and interior-point, both of which are provided by the MATLAB LP toolbox linprog [MATLAB, 2024a]. A detailed comparison of the BD convergence of the four partitioning algorithms using the four datasets is shown in Fig. 15-Fig. 18 in [Qiu, 2024].

As shown in Table 3, both classic LP methods stop prematurely without achieving the optimal solution as they exceed the maximum number of iterations (the parameter MaxIter is set by 10,000). Within the BD framework, on average, km-DV achieves the lowest computation time. This efficiency can be attributed to the balanced sizes of decomposed subproblems, as shown in Table 1 and Table 2. As discussed in Section 3.3, the computation time of BD in each iteration is significantly influenced by the computation time of the largest subproblem, which is relatively smaller when the size variation of decomposed problems is minimized. Additionally, a smaller subproblem, on average, needs fewer iterations to converge to a feasible solution within the BD framework (refer to Fig. 6 in [Qiu, 2024], which demonstrates a positive correlation between the subproblem size and the number of iterations to find a feasible solution).

Finally, we expand the size of the secret dataset from 500 records to 1,000 records and assess the computational efficiency and number of iterations of "k-m-DV" in Fig. 4(a)(b),

• (a) Average computation time
• (b) Avg. number of iterations

Figure 4: Computation efficiency of k-m-DV in larger scale datasets.

respectively. In alignment with the recent related work [Imola et al., 2022], we allocate a calculation time limit of up to 1,800 seconds. Notably, the computation time for "k-m-DV" remains below 6 minutes when using grid maps, text data, and synthetic datasets. However, when using the road map location dataset, the computation time experiences a significant increase, reaching approximately 1,800 seconds when the dataset size is 1,000. This difference can be attributed to the uneven distribution of nodes in the large-scale road network, with low node density in suburban areas and high node density in downtown areas. This leads to some larger decomposed problems, resulting in a longer time to solve.

5 Conclusion

In this paper, we propose improving the scalability of LPbased mDP through secret dataset partitioning and BD. Experimental results with diverse datasets show that our approach expands the scale of existing LP-based mDP by approximately 9 times. Additionally, our findings highlight the importance of balancing the size between decomposed subproblems and MP components when implementing BD.

We identify several promising directions for further improvement. First, while BD enhances metric-DP scalability, its current convergence (up to 650 iterations) is challenging for time-sensitive applications. Leveraging reinforcement learning (RL) could accelerate BD convergence by treating cut selection in Stage 2 as a parameterized stochastic policy. A trained RL model can identify a sequence of cuts for fixed subproblem coefficients, eliminating the need for re-training with each new problem instance. Second, recognizing that when the size of secret datasets is higher (e.g., over 1,000), the decomposed MP might remain too large for current LP solvers to handle. A potential direction is to explore a combination of multiple decomposition techniques, such as DanzigWolfe decomposition, based on the mDP constraint features.

Acknowledgments

This research is partially supported by U.S. NSF grants CNS2136948 and CNS-2313866.