arXiv CS open access

DISTA-LP: Privacy-Preserving Distributed Linear Programming for Metric Differential Privacy

Michael Oluwole et al. · arxiv_cs

arXiv CS · Papers

cryptography cybersecurity privacy security

Linear-programming (LP) mechanisms for metric differential privacy (mDP) enable fine-grained, task-aware utility control, but they are computationally expensive and become privacy-fragile in distributed settings: uploading user-specific LP coefficient matrices (e.g., distance and loss blocks) to a c

Conceptio Open Archive

Analyze & explore related research in Conceptio

Connect this record directly into your research workflow or AI agent. Explore citations across open repositories, query full text, and export structured proof bundles.

Explore in Search →Connect via MCP / API ↗

D ISTA-LP: Privacy-Preserving Distributed Linear Programming for Metric Differential Privacy Michael Oluwole, Ruiyao Liu, Aleksander Stevens, and Chenxi Qiu*

arXiv:2610.27615v1 [cs.CR] 20 Oct 2026

University of North Texas {MichaelOluwole, RuiyaoLiu, AleksanderStevens}@my.unt.edu, Chenxi.Qiu@unt.edu * Corresponding author

utility while maintaining formal guarantees and has proven effective in geo-location [6], [7], [8], [9], natural language [2], [3], and perceptual data such as voice and images [10], [11]. Unlike LDP, designing mechanisms for mDP is complicated by non-uniform privacy constraints and heterogeneous utility sensitivities. Under mDP, the privacy guarantee varies across input pairs, and utility loss depends on both the magnitude and the direction of perturbation [12]. Standard mechanisms, such as the Laplace mechanism [6] and the exponential mechanism [3], satisfy mDP by calibrating noise solely to perturbation magnitude. However, this isotropic approach ignores directional variations in utility across the output space, often resulting in suboptimal utility-privacy trade-offs. To accommodate task-specific utility sensitivities, many recent works turn to optimization-based mechanism design, most commonly formulated as linear programs (LPs) [13], [2], [4]. After discretizing the secret (input) domain X and the reported (output) domain Y, these methods explicitly encode the utility loss of each candidate perturbation and optimize the full perturbation distribution. This yields finely tuned, often near-optimal for the chosen discretization, mechanisms, but shifts the primary challenge to computational scalability: the LP I. I NTRODUCTION optimizes the conditional probabilities P (y | x) for all x ∈ X Local Differential Privacy (LDP) [1] provides a uniform and y ∈ Y, resulting in O(|X ||Y|) decision variables that can indistinguishability guarantee across the entire input domain. quickly reach millions on moderately fine grids, rendering While robust, this one-size-fits-all guarantee can be overly off-the-shelf solvers impractical at scale. conservative for applications that only require local semantic To improve scalability, recent work applies decomposition or spatial obfuscation. For example, in text analytics, enforcing techniques (e.g., Dantzig-Wolfe [12] and Benders [14]) to strong indistinguishability between semantically distant words split a large LP into smaller subproblems, enabling parallelism (e.g., “cancer” vs. “pizza”) can require substantial randomizaand reducing per-solver complexity. Complementarily, hybrid tion, often rendering the released text uninformative; similarly, frameworks [2] integrate the Weighted Exponential Mechanism location privacy may only require hiding a user within a (WEM) with LP by using WEM for inputs where the utilityneighborhood rather than an entire city [2], [3], [4]. loss landscape is relatively flat (i.e., candidate outputs incur To mitigate this utility loss, metric differential privacy similar loss), and reserving LP optimization for utility-critical (mDP) [5] scales indistinguishability via a task-specific metric d. regions where small probability changes can markedly affect By requiring that inputs xi and xj close under d induce similar utility, thereby reducing computation while preserving mDP output distributions, mDP avoids over-perturbation for distant guarantees. Nonetheless, both lines of work still rely on a inputs. This distance-aware framework significantly improves centralized LP backbone, so the maximum feasible secretdomain size |X | remains bounded by single-machine solver capacity (at most ∼ 500 in our setting). More recently, distributed LP frameworks have been proposed to further improve scalability by offloading parts of the Network and Distributed System Security (NDSS) Symposium 2027 optimization to individual users [15], [4]. Specifically, [15] 22–26 March 2027, Seoul, Republic of Korea ISBN 978-1-970672-09-1 lets each user construct a local LP and upload the resulting Abstract—Linear-programming (LP) mechanisms for metric differential privacy (mDP) enable fine-grained, task-aware utility control, but they are computationally expensive and become privacy-fragile in distributed settings: uploading user-specific LP coefficient matrices (e.g., distance and loss blocks) to a central solver can itself leak sensitive information. We present D ISTA-LP, a privacy-preserving distributed mDP optimization framework that replaces raw coefficient uploads with compact surrogate parameterizations and releases only calibrated, perturbed parameter vectors, thereby reducing both upload-stage leakage and problem dimensionality. On the server side, D ISTA-LP reconstructs approximate coefficients and solves a scalable mastersubproblem formulation via Benders decomposition. To certify end-to-end privacy, we introduce a probabilistic mDP (PmDP) auditor that empirically quantifies the privacy loss induced by parameter release, allocates the required privacy budget prior to optimization, and composes it with the downstream mechanism to obtain an end-to-end mDP guarantee. Experiments on large realworld road-network datasets (Rome, New York City, and London) show that D ISTA-LP achieves near-zero mDP violation rates, substantially reduces computation time compared to monolithic and hybrid baselines, and matches or improves utility, thereby demonstrating stronger upload privacy and significantly improved scalability.

https://dx.doi.org/10.14722/ndss.2027.230054 www.ndss-symposium.org

coefficient blocks, after which the server solves a global LP that aggregates all users’ contributions. In parallel, [4] proposes an anchor-based two-phase framework in which each user selects a small set of anchor records, allowing the server to solve a compact LP over a reduced domain. While these designs substantially reduce runtime and retain competitive utility, their privacy treatment remains limited: [15] does not account for potential privacy leakage from the uploaded coefficient matrices, and [4] must spend additional privacy budget to compensate for anchor-induced distance approximation errors, which can increase utility loss.

best fidelity–privacy trade-off. Moreover, we design the parameter perturbation to be structure-preserving and to minimize surrogate reconstruction error, so that the induced coefficient approximation does not degrade utility or violate LP feasibility. Because the uploaded parameters lie in a continuous space, tightly bounding the induced privacy-loss distribution in closed form can be challenging under structure-preserving perturbations. We therefore adopt probabilistic mDP (Definition 3) accounting and certify upload-stage privacy using a samplingbased posterior leakage auditor (based on Proposition 2). The server then reconstructs approximate coefficient blocks from the perturbed parameters and solves the subsequent LP using the Benders decomposition framework.

Our Contributions

Contribution 1: Distributed Approximation LP (D ISTA-LP) Contribution 3: Empirical validation on large-scale geframework. In this paper, we propose D ISTA-LP, a distributed olocation datasets. We conduct extensive experiments on LP framework that optimizes the perturbation distribution large-scale real-world road-network datasets from Rome, New and enables scalable perturbation optimization under mDP. York City, and London [20] to validate the effectiveness of In D ISTA-LP, each user m locally constructs a set of LP D ISTA-LP. Our results show that D ISTA-LP substantially coefficient blocks that capture the distance and utility-loss improves scalability over centralized LP-based mechanisms, structure within its relevance neighborhood Nm (Nm ⊂ X reducing solver runtime and memory consumption while and |Nm | ≪ |X |). A key distinguishing feature is that we maintaining competitive, and often superior, utility. Compared explicitly model and quantify the privacy loss incurred by with predefined-noise baselines (Laplace [6] and the exponential coefficient uploads, an aspect overlooked in prior distributed mechanism [21]), LP-based methods [13], [4], and hybrid mDP formulations such as [15]. Because these coefficient approaches such as COPT [2], D ISTA-LP achieves lower blocks are derived from users’ authentic records, they can expected utility loss while attaining near-zero empirical mDP inadvertently reveal sensitive spatial or semantic correlations. violation ratios under comparable privacy budgets. We further To mitigate this risk, D ISTA-LP analytically characterizes the demonstrate that the proposed surrogate-based coefficient apprivacy loss induced by coefficient uploads and allocates a proximation enables communication- and computation-efficient corresponding portion of the global privacy budget using the uploads, and that the probabilistic auditing procedure provides sequential composition property of mDP [4], ensuring that the accurate empirical estimates of upload-stage privacy leakage. Collectively, these findings indicate that D ISTA-LP makes cumulative leakage remains rigorously bounded. On the server side, D ISTA-LP uses the uploaded LP LP-based mDP optimization practical at city scale without coefficients to assemble a global LP that integrates a WEM sacrificing formal privacy guarantees. Section II reviews the fundamentals of mDP, defines privacy structure [2], [15], handling distant record pairs analytically loss, and introduces the baseline LP formulation. Section III while explicitly optimizing over nearby, utility-sensitive ones. presents an overview of our approach and highlights the key reTo further scale, D ISTA-LP adopts Benders decomposition: search challenges. Section IV describes the distributed pipeline, a compact master problem optimizes the global mechanism including surrogate construction, parameter perturbation, and variables and WEM weights, while each user solves an reconstruction, along with the WEM+LP formulation enhanced independent neighborhood subproblem over Nm that acts as a by Benders decomposition. Section V reports the experimental feasibility/optimality check and returns cutting planes (from evaluation of our approach. Section VI discusses related work, the subproblem dual) to iteratively tighten the master [14]. and Section VII concludes the paper. Contribution 2: Coefficient surrogate modeling, parameter perturbation, and probabilistic mDP accounting. Notably, II. P RELIMINARIES LP coefficient blocks are deterministically derived from users’ In this section, we present the preliminaries of mDP, true records; thus, uploading them directly can leak inforincluding the underlying threat model, the formal definition of mation about the underlying records. To mitigate this risk, mDP, the LP-based optimization approaches, and the WEM we replace raw coefficient uploads with a surrogate-based integration. approach: each user fits a low-dimensional surrogate that approximates the required coefficient blocks and uploads only Threat model. In general, a data perturbation mechanism can perturbed surrogate parameters. We consider both low-rank be represented as a probabilistic mapping M : X → Y, where models (truncated singular value decomposition (SVD)) and the domain X denotes the users’ secret data domain, and the parametric/basis-function surrogates (polynomial, Gaussian, range Y denotes the perturbed data domain. Based on the and radial basis function (RBF)) [16], [17], [18], [19], and intrinsic characteristics of X , a metric d : X 2 → R is defined systematically compare them in terms of approximation fidelity to quantify the distance between any two records xi , xj ∈ X , and privacy leakage; in our setting, low-rank SVD offers the denoted as dxi ,xj .

We adopt a standard honest-but-curious threat model [22], where the server faithfully executes the protocol but attempts to infer sensitive information from received reports. Specifically, given the perturbation mechanism M and a prior distribution of X, the server can apply Bayes’ rule to compute the posterior distribution of X upon observing a perturbed record y [7]. This model reflects realistic deployments where datacollecting organizations are not fully trusted. Conversely, we assume users are honest and correctly apply the perturbation protocol; collusion, either among users or between users and the server, is excluded as it falls outside the scope of LDP/mDP frameworks [23]. To formally bound the server’s ability to distinguish between secrets via these posterior distributions, we define the privacy loss as follows:

Column 𝑘

Follow linear mDP constraints

Perturbation probability

Follow exponentialdecay form

1 Pr[M(xi ) = y] . ln dxi ,xj Pr[M(xj ) = y]

Distant, low-impact utility Distance region to 𝑦𝑘

Fig. 1: Idea of WEM+LP. by a stochastic matrix Z = {zi,k }(xi ,yk )∈X ×Y , where zi,k denotes the probability of outputting yk given input xi . Under this discretization, the mDP condition in Eq. (2) can be written as the following linear constraints:

(1)

zi,k − eϵdxi ,xj zj,k ≤ 0,

Intuitively, PLxi ,xj (y) quantifies how differently the mechanism can produce output y when the input is xi rather than xj , normalized by the metric distance between xi and xj . Smaller values indicate stronger indistinguishability between xi and xj , and hence less information leakage to an adversary (e.g., the server) through posterior inference.

∀xi , xj ∈ X , ∀yk ∈ Y.

(3)

Each row of Z must also satisfy the normalization constraints: PK ∀xi ∈ X . (4) k=1 zi,k = 1, Let cxi ,yk denote the utility loss incurred when reporting yk given the true record xi , and let πxi denote the prior probability that the true record equals xi . Define the cost matrix Acy ∈ R|X |×|Y| by [Acy ]i,k := πxi cxi ,yk . Then the expected utility loss is PN PK L(Z; Acy ) = i=1 k=1 [Acy ]i,k zi,k . (5)

In the remainder of the paper, for brevity, we use “privacy loss” to refer to “metric-normalized privacy loss.” Definition 2 (mDP [6]). We say that a perturbation mechanism M : X → Y satisfies (ϵ, d)-mDP if and only if that the privacy loss is uniformly bounded by ϵ, i.e., sup sup PLxi ,xj (y) ≤ ϵ.

ℐ(𝑦𝑘 )

…

ℐc (𝑦𝑘 )

ℐ(𝑦𝑘 )

Definition 1 (Metric-normalized privacy loss). Let (X , d) be a metric space of secrets, and let M be a perturbation mechanism with input domain X and output space Y. We define the metricnormalized privacy loss between two inputs xi , xj ∈ X at output y ∈ Y as PLxi ,xj (y) =

Near, highimpact utility regions

…

The optimal perturbation matrix is obtained by the LP min L(Z; Acy )

(2)

xi ̸=xj y∈Y

s.t. Eqs. (3) and (4) hold.

(6)

This formulation involves O(|X ||Y|) decision variables and O(|X |2 |Y|) constraints, and the pairwise mDP constraints quickly become the bottleneck for efficiency.

Over more than a decade of research on mDP, a variety of mechanisms have been proposed to satisfy mDP, including Laplace-based approaches [6] and exponential mechanisms [24]. While these predefined mechanisms offer general-purpose solutions with clear theoretical guarantees, their design is largely driven by distance alone, and therefore they generally do not capture the heterogeneous sensitivity of utility loss under perturbations in different directions. As a result, they may be suboptimal in terms of privacy-utility tradeoff. This observation has motivated optimization-based approaches, which seek to directly construct perturbation mechanisms that minimize utility loss subject to mDP constraints.

Weighted exponential mechanism (WEM) integration. To improve scalability, prior works such as [2], [15] impose a WEM-inspired structure on many perturbation probabilities. As illustrated in Fig. 1, for each fixed output yk , the input domain X is partitioned into a nearby, high-impact utility region I(yk ) and a distant, low-impact utility region I c (yk ). Instead of optimizing every probability zi,k through LP, the probabilities corresponding to xi ∈ I c (yk ) are constrained to share an exponential-decay form. The rationale is that distant input-output pairs typically reside in relatively flat regions Optimization (LP)-based approaches. Optimizing a per- of the utility landscape and thus have limited utility impact, turbation mechanism M over continuous secret or output making a structured exponential tail sufficient. spaces is computationally challenging. Therefore, a common Concretely, let w = [w1 , . . . , wK ] be nonnegative output approach is to discretize X and Y into finite sets [2]. Let weights. For each output yk , define an exempt set I(yk ) ⊆ X , K X = {xi }N for example, the r nearest xi to yk [2] or those within a distance i=1 and Y = {yk }k=1 . Then, M can be represented

0.06

0.014

Privacy budget = 4.0

Feasible region: 𝛀 Server

dy dx dy 𝐀dx 1 … , 𝐀 𝑀 , 𝐀1 … , 𝐀 𝑀

Upload

Probability mass

Objective: Minimize ℒ 𝐙; 𝐀1 … , 𝐀𝑀

0.01 0.008 0.006 0.004

Privacy loss Privacy loss >

0.002

Coefficients dy cy 𝐀dx 1 , 𝐀1 , 𝐀1

User 1

…

Coefficients dy cy 𝐀dx 𝑀 , 𝐀𝑀 , 𝐀𝑀

0.04 0.03 0.02

100

150

0.02

0.015

0.01

0.005

Privacy loss Privacy loss > 60

Privacy loss Privacy loss >

0 0

200

Privacy loss

(a) Adx

User M

Fig. 2: Distributed LP Framework: Uploaddy cy ing Coefficient Blocks Adx m , Am , Am by each user m (m = 1, ..., M ).

0.05

0.01

0 0

Privacy budget = 4.0

0.012

Probability mass

Linear programming

400

600

800

Privacy loss

(b) Ady

dy Fig. 3: Empirical distributions of the privacy loss induced by uploading Adx m , Am , cy and Am . For ϵ = 4, the fractions of samples exceeding the budget are 96.97%, 80.64%, and 79.73%, respectively (exceedances highlighted in red).

threshold [15]. For every non-exempt (“far”) pair xi ∈ / I(yk ), we impose the WEM constraints zi,k = wk exp − 2ϵ d(xi , yk ) , ∀yk ∈ Y, ∀xi ∈ / I(yk ). (7)

where (Z, w) is obtained by solving FLP+WEM Acy , Adx , Ady =

Hence, the optimization no longer treats all |X | |Y| probabilities as independent PK decision variables; only the exempt probabilities, totaling k=1 |I(yk )|, together with the K weight variables in w, remain to be optimized, while the row-normalization constraints in Eq. (4) distribute the remaining probability mass among the exempt outputs. To represent the WEM-integrated formulation compactly, we collect three coefficient blocks. Let T := {dx, dy, cy} index these blocks, and denote them collectively by {At }t∈T , where

arg

min

(Z,w)∈Ω(Adx ,Ady )

L (Z; A ) .

(13) (14)

III. F RAMEWORK M OTIVATION AND C HALLENGES

Notably, although WEM constraints reduce the effective search space by enforcing an exponential-decay form for many “far” transition probabilities, its computational cost remains prohibitive on large-scale domains. This blow-up is largely driven by the mDP requirement, which introduces distance-based indistinguishability constraints for (potentially) every pair of secrets in X . Although these constraints are dx |X |×|X | dx A ∈R , [A ]i,j := dxi ,xj , (8) defined globally, their effect is inherently local: as dxi ,xj increases, the corresponding constraint becomes increasingly Ady ∈ R|X |×|Y| , [Ady ]i,k := dxi ,yk , (9) permissive, so distant records exert much less influence than Acy ∈ R|X |×|Y| , [Acy ]i,k := πxi cxi ,yk . (10) nearby neighbors [2]. This locality motivates scalable LPHere, dx captures pairwise distances within X and determines based methods that focus on enforcing indistinguishability how strongly different input rows are coupled by the mDP within local neighborhoods, rather than explicitly instantiating constraints. The block dy records input-output distances and is constraints for all secret pairs. Hence, a promising direction is used to define the near/far structure around each output yk as to apply a distributed framework that offloads local structure well as the exponential-decay terms in the WEM constraints. to users. Fig. 2 illustrates a distributed LP framework adopted in Finally, cy encodes the prior-weighted utility loss of reporting prior work (e.g., [15]). Each user m constructs local coefficient yk when the true input is xi , and therefore supplies the dy cy blocks Adx m , Am , Am restricted to a neighborhood Nm around coefficients of the utility objective. With WEM constraints, the feasible region is the convex the true record xm (e.g., the 20 nearest secrets to xm ), rather than spanning the entire secret domain X . The server aggregates polyhedron these blocks across users and solves the resulting integrated Ω Adx , Ady (11) LP of the following form (the full formulation is deferred to   dx Appendix C for completeness): zi,k − eϵ[A ]i,j zj,k ≤ 0,     2   ∀(x P i , xj , yk ) ∈ X × Y, (Z, w) min L(Z; Acy , ..., Acy (15) M) = . |X |×|Y| × RK 1 yk ∈Y zi,k = 1, ∀xi ∈ X , ∈ [0, 1] dy 1 ≥0   dy dy dx dx  s.t. Z ∈ Ω A1 , ..., AM , A1 , ..., AM . (16) zi,k − wk e− 2 ϵ[A ]i,k = 0,    ∀yk ∈ Y, ∀xi ∈ I c (yk )

This integrated program retains the constraint template of the centralized LP FLP+WEM ((At )t∈T ) as described in Eq. (13), but restricts it to pairs within each local neighborhood Nm , resulting in a much smaller optimization problem. Accordingly, the number of decision variables drops from O(|X ||Y|) to

LP+WEM generator as an operator. We view the LP-based mechanism generator with WEM integration (LP+WEM) as a functional operator FLP+WEM : At t∈T 7→ (Z, w) ∈ [0, 1]|X |×|Y| × RK ≥0 , (12)

dy cy , and the number of constraints from Approximating Adx m , Am , Am with low-dimensional surrogates P inevitably introduces reconstruction error, which can affect M 2 O(|X |2 |Y|) to O m=1 |Nm | |Y| , providing substantial both the feasibility of mDP constraints and the utility of PM scalability gains when m=1 |Nm | ≪ |X |. the LP solution produced by FLP . In particular, distance However, the distributed design in [15] implicitly treats overestimation can inflate terms such as eϵd , enlarge the feasible dy cy the upload of local LP coefficient blocks Adx region, and yield mechanisms that violate the intended mDP m , Am , Am as privacy-neutral. This assumption is problematic because these constraints. Therefore, we must design calibrated asymmetric blocks are deterministic functions of users’ authentic records perturbation and conservative reconstruction procedures so and can themselves reveal sensitive information about the that the recovered coefficients remain mDP-preserving while underlying data. To formalize this, we define the coefficient- retaining high utility. upload mechanism Mtc as the mapping from a true record Challenge 3: How can we evaluate and account for the x ∈ X to the uploaded coefficient block for type t ∈ T : privacy loss of releasing continuous surrogate parameters? e t lie in continuous spaces. Mtc (xm ) ≜ Atm . (17) The released surrogate parameters Θ m For the corresponding parameter-noising mechanisms, deriving t If no randomization is applied, then Mc is deterministic. closed-form privacy-loss bounds or privacy-loss distributions is Consequently, for two distinct records xi ̸= xj , there exists an generally difficult, especially under asymmetric, truncated, and e (e.g., A e = At )) such that uploaded observation A i otherwise nonstandard perturbations. Hence, we need practical h i h i t t evaluation and accounting methods to quantify the induced e e Pr Mc (xi ) = A = 1 while Pr Mc (xj ) = A = 0. (18) privacy loss and incorporate it into end-to-end budget allocation and composition. e The corresponding privacy loss for observing A, Next, we present our method and explain how it addresses ! t e each of the three challenges. (x ) = A] 1 Pr[M c i e = , (19) PLxi ,xj A ln t e dxi ,xj Pr[Mc (xj ) = A] IV. M ETHODOLOGY is infinite, showing that uploading coefficient blocks without In this section, we detail the design of DistA-LP and randomization can cause unbounded leakage and thus requires explain how it addresses Challenges 1–3 in Section III. As explicit protection at the upload stage. Fig. 4 shows, DistA-LP decomposes end-to-end mechanism Notably, simply adding naive noise to the coefficient blocks generation into three coordinated steps: Atm (for t ∈ {dx, dy, cy}) remains insufficient. As a base- ⃝ 1 Coefficient surrogatization. Each user m locally conline, we perturb each block entrywise with Laplace noise structs three coefficient blocks: the intra-neighborhood Lap 0, ∆ϵt (with ϵ = 4), where ∆t is the metric sensitivity distance matrix Adx m , the neighborhood-output distance of Atm under d(·, ·) (see Appendix C). We then measure the matrix Ady , and the utility-loss matrix Acy m m . To reduce e induced privacy loss PLxi ,xj A . Figure 3(a)–(c) shows that information leakage, user m fits a compact coefficient sura substantial fraction of the mass still lies above the target rogate f (·; Θtm ) (parameterized by Θtm ) to approximate budget. Atm , i.e., [Atm ]i,j ≈ f (i, j; Θtm ), t ∈ {dx, dy, cy}. Motivated by these examples, we propose the D ISTA-LP 2 Private parameter release. The user privatizes the ⃝ framework to explicitly account for the privacy loss incurred by surrogate parameters using a local randomizer Mpt , coefficient uploads and to allocate a dedicated budget for correct e t = Mp (Θt ), and uploads Θ e t (t ∈ {dx, dy, cy}) Θ t m m m end-to-end composition, raising several unique challenges to to the server. address: 3 Perturbation Optimization. The server reconstructs each ⃝ Challenge 1: How can we reduce the privacy loss incurred e t from the noised surapproximate coefficient block A m by uploading local LP coefficient blocks? As illustrated e t ]i,j = f (i, j; Θ e t ), aggregates rogate parameters via [ A m m in Fig. 3, directly uploading coefficient blocks (even with the reconstructed blocks, and solves the resulting LP to naive perturbation) can leak highly sensitive information ∗ obtain each optimal local perturbation matrix Z . N about a user’s locally relevant neighborhood Nm and, by m ×Y t Finally, each user downloads the resulting perturbation matrix extension, the underlying record. A key reason is that Am is high-dimensional and encodes fine-grained, instance-specific and applies it to map its true record xm to a randomized output neighborhood structure. To mitigate this leakage, we replace ym . The server then uses the perturbed records for downstream each high-dimensional block Atm with a low-dimensional analytical tasks. surrogate representation and then apply a parameter-noising Privacy accounting across uploads. As shown in Fig. 4, e t , aiming to suppress D ISTA-LP involves four user-to-server releases: three uploads mechanism to the surrogate parameters Θ m fine-grained, instance-specific structure and reduce record-level of the privatized surrogate parameters Θ e dx , Θ e dy , and Θ e cy in m m m distinguishability. 2 followed by a final upload of the perturbed report ym Step ⃝, Challenge 2: How can we minimize the impact of surro- after applying the downloaded perturbation matrix. Since each gate approximation on both mDP feasibility and utility? release consumes privacy budget, we assign per-release budgets

M m=1 |Nm ||Y|

Server

Optimal perturbation matrix

Step ③ Perturbation optimization ෩cy ෩cy Objective function: 𝐀 1 , … , 𝐀𝑀 dy

෩dx ෩dx ෩ ෩ mDP constraints 𝐀 1 , … , 𝐀 𝑀 , 𝐀1 , … , 𝐀 𝑀 Reconstructed coefficients

𝜺𝐝𝐱

෩cy 𝐀 𝑚

෩dy 𝐀 𝑚

෩dx 𝐀 𝑚 Upload

Privacy budget:

solution

𝜺𝐝𝐲

𝜺𝐜𝐲

Step ② Private parameter release

෩ dx 𝚯 𝑚

෩ dy 𝚯 𝑚

෩ cy 𝚯 𝑚

Step ① Coefficient surrogatization

𝚯dx 𝑚

𝚯𝑚

𝐀dx 𝑚

𝐀𝑚

cy 𝐀𝑚

User 𝑚

𝒁𝒩𝑚 ,𝒴

Upload

Downstream task Upload

≤ 𝜺

𝜺𝐩𝐭

𝑦𝑚

Download

Record perturbation

Real record 𝑥𝑚 Fig. 4: The DistA-LP framework.

e dx , Θ e dy , Θ e cy , and ym , respectively, ϵdx , ϵdy , ϵcy , and ϵpt to Θ m m m and enforce the end-to-end constraint ϵdx + ϵdy + ϵcy + ϵpt ≤ ϵ,

෡𝑡𝑚 Permutated coefficient 𝐀

Approximated surrogate parameters 𝚯𝑡𝑚

(20)

where ϵ is the predetermined overall privacy threshold. Since 3 are the server-side reconstruction and LP solving in Step ⃝ pure post-processing of the uploaded privatized parameters, they incur no additional privacy cost; therefore, tracking and composing these four budgets suffices to guarantee the end-toend privacy level of the full protocol. 1 ⃝ 3 in Sections IV-A– Next, we present the details of Steps ⃝– IV-C, respectively.

① Approximation

② Perturbation

1 Coefficient Surrogatization A. Step ⃝: ③ Reconstruction

In this step, we approximate each coefficient block Atm (t ∈ {dx, dy, cy}) using a compact parametric surrogate f (·; Θtm ), instantiated either as a predefined functional form or a lowrank SVD approximation, with parameters Θtm . Restricting Atm to a low-dimensional model class suppresses fine-grained, instance-specific variations within Nm that could otherwise be revealing. In addition, the mapping Atm 7→ Θtm is typically many-to-one, so multiple distinct blocks can induce the same (or similar) parameter vector; as a result, Θtm is generally less distinguishable than Atm itself, increasing an adversary’s uncertainty about the underlying coefficients. We next detail the two approximation families: predefined parametric functions and SVD-based surrogates. Local permutation. To reduce the similarity among coefficient blocks across users, we apply a local permutation to the ordering of locations within each neighborhood Nm . Let Em ∈ {0, 1}|Nm |×|Nm | be the corresponding permutation matrix, with exactly one 1 in each row and each ⊤ column (and 0 elsewhere); hence E⊤ m Em = Em Em = I −1 ⊤ and Em = Em . The permutation for each user can be

෩𝑡𝑚 Reconstructed coefficient 𝐀

Perturbed surrogate ෩ 𝑡𝑚 parameters 𝚯

Fig. 5: Surrogate-based coefficient privatization pipeline (permutation, approximation, parameter perturbation, and reconstrucet → A et . tion): Atm → Θtm → Θ m m

generated via a deterministic procedure (e.g., the optimization in Eq. (17)), so the coordinator may know (or be able to reproduce) the set {Em }M m=1 . However, surrogate/coefficient uploads are received without user identifiers in the anonymized upload phase; thus, upon receiving a permuted block, the coordinator cannot determine which Em was applied and cannot undo the permutation. Notably, Em permutes only the neighborhood index in Nm

TABLE I. Comparison of different matrix approximation methods. Surrogate 2D Gaussian

Size of Θtm 6

Accuracy Low

2D Polynomial

Medium

RBF

High

Low-Rank SVD

30–50

Very High

Notes Compact, but only fits symmetric bell-shaped data well. Simple to fit; more effective for smooth surface trends. Flexible and localized, better for multi-peak or irregular matrices. Captures global structure with fewer components; efficient to store and reconstruct.

and leaves the output index set Y unchanged. Given Em , we permute each coefficient block Atm as ( t ⊤ t b = Em Am Em , t = dx, A (21) m t Em A m , t ∈ {dy, cy},

Surrogate f (·; Θtm ) derivation. Each user then approximates b t (t ∈ {dx, dy, cy}) using the their permuted blocks A m parametric surrogate f (·; Θtm ): ∀m ∈ [M ], ∀t ∈ T .

(22)

Here, we instantiate f (·; Θtm ) using four parametric families that trade expressiveness for parameter count (Table VII). We jointly determine the permutations {Em }M m=1 and surrogate parameters {Θtm }m∈[M ], t∈T by minimizing the total entrywise ℓ1 fitting error: min

{Em ,Θtm }M m=1

M X X m=1 t∈T

b t − f (·; Θt ) λt A m m

, 1

Subtract noise from singular values (Sk ) or from the reconstructed matrix directly to ensure conservative output.

entrywise ℓ1 reconstruction loss in (23) using the Nelder– Mead simplex algorithm [25], a robust, derivative-free optimizer well-suited for nonconvex objectives and nonsmooth losses. (2) Permutation search. With {Θtm }M m=1 fixed, we refine the local permutation Em via greedy pairwise index swaps (i, j) that decrease the global objective in (23). We accept only swaps that reduce the total weighted ℓ1 loss, ensuring monotonic descent. After each accepted swap, we refit {Θtm }M m=1 under the updated ordering for consistency. According to the experimental results in Fig. 6–7 this alternating scheme converges to a locally optimal configuration in practice and typically stabilizes within 10-20 iterations. Table VII compares surrogate parametric families in terms of parameter complexity, approximation accuracy, and the noise strategies used to ensure conservative, privacy-preserving outputs. Overall, 2D Gaussian and polynomial surrogates offer compact representations but are best suited to smooth, globally structured surfaces (and, in the Gaussian case, near-symmetric shapes), whereas RBFs provide greater expressiveness at the cost of substantially larger parameter sets. In contrast, the low-rank SVD surrogate provides the most favorable balance, preserving global structure with high fidelity while keeping the representation low-dimensional.

i.e., for the square block (t = dx) indexed by Nm × Nm , we permute both neighborhood axes, while for the rectangular blocks (t ∈ {dy, cy}) indexed by Nm × Y, we permute only the neighborhood axis. We do not permute Y because the server needs to formulate a single optimization problem in a common output index space, so that decision variables associated with y ∈ Y are aligned across users; otherwise, cross-user consistency of output-indexed 3 in Section IV-C). variables would be ill-defined (see Step ⃝ After user m downloads the optimized local perturbation matrix (indexed by the permuted neighborhood order), it applies the inverse permutation locally to recover the original indexing ⊤ (i.e., left-multiplies by E−1 m = Em ).

b t ≈ f (·; Θt ), A m m

Noise Strategy Subtract Laplace noise from amplitude (A) and increase spread (σ) to lower the surface everywhere. Subtract Laplace noise from coefficients, especially the intercept, to bias the surface downward. Add negative-biased noise to weights (αi − |Lap(b)|); optionally perturb centers for robustness.

2 Private Parameter Release B. Step ⃝: To preserve privacy in our distributed LP, we introduce three randomized parameter-perturbation mechanisms Mpt (t ∈ {dx, dy, cy}) that inject calibrated and asymmetric noise into the uploaded surrogate parameters: e t = Mp Θt , Θ ∀m ∈ [M ], ∀t ∈ T . (24) t m m

(23)

where λt ≥ 0 controls the trade-off among approximation components, and we fix λdx = 1 without loss of generality. Here, we prefer ℓ1 over ℓ2 for robustness to outliers and discretization noise, which typically yields more stable fitting.

Here, each Mpt is a randomized mapping from the surrogateparameter space to itself,

Joint optimization of permutation {Em , Θtm }M m=1 . We adopt

where dt is the number of parameters in Θtm for block type

an alternating scheme between (1) fitting parameters for a fixed permutation and (2) improving the permutation through greedy swaps: (1) Parameter fitting. Holding {Em }M m=1 fixed, for each surrogate family we fit Θtm by minimizing the weighted

t. This is analogous to the perturbation mechanisms defined earlier for raw records, but with a different domain: instead of mapping secrets x ∈ X to perturbed outputs y ∈ Y, Mpt maps e t prior to surrogate coefficients Θtm to noised coefficients Θ m upload.

Mpt : Rdt → Rdt ,

(25)

(1)(2) Mechanisms Mpdx and Mpdy . Mechanism Mpdx perturbs the surrogate parameters Θdx m for the intra-neighborhood distance block. Mechanism Mpdy perturbs Θdy m for the crossdomain distance block, which links each local neighborhood to the shared output space Y. Both Mpdx and Mpdy add calibrated Laplace noise with an asymmetric form, e dx = Θdx − |η|, Θ m m

the privacy-loss bound is violated. Equivalently, the bound PLxi ,xj (z) ≤ ϵ holds with probability at least 1 − δ for a uniformly selected input pair and its corresponding randomized output. When δ = 0, the definition recovers exact PL as given in Eq. (2). Proposition 1. (Sequential composition of PmDP, adapted from [4]) Let {Mpt }t∈T be randomized mechanisms with conditionally independent outputs P given the input. If each Mpt satisfies (ϵt , δt )-PmDP and t∈T δtP < 1, thenPthe joint mechanism M = (Mpt )t∈T satisfies ( t∈T ϵt , t∈T δt )PmDP.

η ∼ Lap(0, b),

and truncate the perturbed values to a bounded range to ensure numerical stability and bounded sensitivity. Notably, the onesided noise is intended to bias the perturbation toward reducing the distance values, thereby tightening (rather than relaxing) the corresponding mDP constraints and avoiding infeasible solutions that could cause mDP violations.

Proposition 2. [Concentration Guarantee for Probabilistic mDP Sampling] If the empirical violation hrate satisfies i p̂Sℓ = ξδ for some constant ξ < 1, then Pr pXℓ2 ≤ δ ≥ 1 − 2 exp(−2Sℓ (1 − ξ)2 δ 2 ). The detailed proof can be found in Appendix B.

(3) Mechanism Mpcy (utility-loss parameters). Mpcy perturbs the surrogate parameters for the utility-loss block Θcy m used in the LP objective. This prevents the released objective coefficients from inadvertently revealing sensitive information encoded in the utility-loss structure (e.g., user-specific preferences or constraints), while keeping the perturbed coefficients well-behaved so that the global LP remains feasible and numerically stable in practice, we use c̃i,k = ci,k + |ηc |,

3 Perturbation Optimization C. Step ⃝:

ηc ∼ Lap (0, bc ) .

The detailed description of perturbation mechanisms for both predefined parametric surrogates and SVD-based approximations can be found in Appendix C. a) Privacy Loss Evaluation.: According to Definition 1, t given the observation of perturbed coefficient surrogate Θ̃m , the privacy loss (PL) of any two records xi , xj ∈ X is defined i by  h p t Pr M (x ) = Θ̃ i t m 1 t i  . (26) PLxi ,xj (Θ̃m ) = ln  h t p dxi ,xj Pr M (x ) = Θ̃ t

The server-side after perturbation, depending on the surrogate family used during fitting, the server reconstructs each coefficient block through a predefined parametric mapping or through a low-rank factorization and then applies a structurepreserving post-processing step (which is privacy-free due to the post-processing immunity) to ensure feasibility and improve downstream stability. Server-side reconstruction (parametric and low-rank surrogates). The coordinator receives the privatized surrogate e t }t∈T from each user m and reconstructs each parameters {Θ m coefficient block in the permuted index space used during local fitting. The reconstruction follows the corresponding surrogate t t e e family: Am = f ·; Θm , t ∈ T .

Note that Θ̃m is defined over a continuous domain. Because we deliberately use non-standard (asymmetric and truncated) parameter noise to reduce the violation rate, the distribution t of PLxi ,xj (Θ̃m ) does not admit a closed-form expression. For computational tractability, we therefore estimate the privacy loss via Monte Carlo sampling. Here, we define pX 2 as the probability PL is violated: pX 2 = Pr PLXi ,Xj (Z) > ϵ , (27)

∁ For each user m, we define Hm and Hm to denote the sets of real-perturbed location pairs that follow and do not follow the WEM constraints, respectively. Specifically,

where Xi , Xj are random variables representing two secret records in X , we require pX 2 ≤ δ, i.e., the posterior leakage bound can be satisfied with a high probability 1 − δ.

(xi , yk ) ∈ Nm × Y | zi,k − wk e−ϵdxi ,yk = 0

∁ Hm

(Nm × Y) \ Hm .

Definition 3. (ϵ, δ)-Probabilistic mDP (PmDP). A randomized mechanism Mℓ : X → Z satisfies (ϵ, δ)-PmDP if X 1 Pr PLxi ,xj (z) > ϵ ≤ δ, (28) 2 |X | z∼Mℓ (xi )

Moreover, given each perturbed record yk ∈ Y, we define Hm,k = {xim ∈ Nm | (xim , yk ) ∈ Hm } as the set of user m’s peer records selecting yk by following WEM.

where δ ∈ [0, 1) is the average probability, over uniformly selected input pairs (xi , xj ) and the randomness of Mℓ , that

Following prior work [2], we apply a heuristic strategy to determine Hm given User m’s uploaded distance matrix D̃Nm ×Y and utility loss matrix C̃Nm ×Y .

xi ,xj ∈X

constraints: CHm zHm ≥ bNm − CH∁m z∗H∁ . If the guessed m value is suboptimal or infeasible, the subproblem generates a new constraint (or cut) to refine the feasible region V, improving the MP’s future estimates of v. The constraint matrices CHm across different subproblems are disjoint, enabling parallel validation of the guessed values of vm . This iterative process continues, with each subproblem progressively narrowing the solution space until convergence to the optimal solution is achieved. Finally, the coordinator restores the original ordering by applying the inverse permutation. For the square intraneighborhood distance block,

We rewrite the objective function as M X

L (ZNm ) m=1 PM P m=1 (xi ,yk )∈Hm ci,k zi,k |

m=1

}

ci,k zi,k (xi ,yk )∈H∁ m

| =

Following WEM

Not following WEM

k=1 αk wk +

(29)

(30)

}

z ∁ m=1 cH∁ m Hm

(31)

ϵpt d(xi ,yℓ ) PM P 2 where each αk = m=1 xi ∈Hm,k ci,l e− is a constant. The constraints described in Eq. (3)(4) and Eq. (7) can be rewritten as

CH zHm | m {z }

Following WEM

CH∁ zH∁m {z } | m

≥ bNm

b m. Zm = Em Z

(32)

(36)

D. A Two-User Running Example

Consider two users with local neighborhoods N1 = {x1 , x2 , x3 } and N2 = {x4 , x5 , x6 }. and a shared public output where CHm , CH∁m includes the mDP constraints between domain Y = {y1 , y2 , y3 }. Suppose User 1’s true location is x2 the perturbation vectors of the records in Hm : and User 2’s is x5 . Neither true location is sent to the server h i during mechanism generation. CHm , CH∁m Local construction and release. Each user m constructs three   blocks: Adx .. . m , containing distances between locations in Nm ;  . .. dy · · · · · · · · · A , containing distances from Nm to Y; and Acy   m m , containing ∀(x , y), ϵpt dxi ,xj i  ···  1 · · · −e · · · the corresponding utility losses. All three blocks are 3 × 3 in  (x , y) =  j  · · · −eϵpt dxi ,xj · · · 1 ···    ∈ Nm this example. The user permutes the neighborhood indices, fits  . .. a surrogate to each permuted block (for example, a rank-r SVD .. . ··· ··· ··· approximation), perturbs the fitted parameters, and uploads the Here, the privacy budget assigned to the perturbation matrix three perturbed parameter sets under budgets ϵdx , ϵdy , and ϵcy . optimization is ϵpt . By sequential composition for PmDP The user retains the mapping between the permuted rows and (Proposition 1), if ϵpt ≤ ϵ − (ϵdx + ϵdy + ϵcy ), then the overall the original locations. mechanism satisfies the target (ϵ, δ)-PmDP guarantee. Reconstruction and accounting. From each upload, the server Like [15], we apply Benders’ decomposition to decompose reconstructs an approximate 3 × 3 block. The privacy auditor the LP into two simpler stages: a master program and a set of evaluates the disclosure associated with the parameter releases subproblems. and determines whether their allocated budgets meet the stated Stage 1 - Master Program (MP): The MP focuses on privacy criterion. The remaining budget ϵpt is used for the determining the weights of EM w = [w1 , . . . , wK ]. To simplify output mechanism, with equation 20. Reconstruction and the optimization, it replaces the utility loss terms associated optimization operate on the released parameters and require with data perturbation for each user m, given by cH∁m zH∁m , no further upload of the raw blocks. with a single decision variable vm , i.e., vm = cH∁m zH∁m . The Optimization. Using the reconstructed distance and loss blocks, MP is formulated as the following LP problem: the server identifies the entries handled by the WEM structure PK and the entries optimized in the local LPs. Its Benders master min αk wk (33) problem chooses the shared WEM weights and an estimated | k=1{z } Utility loss (following WEM) cost for each user. Subproblem 1 checks the master solution PM against the constraints for N1 ; Subproblem 2 does the same for + v (34) m | m=1 {z } N2 . An infeasible or underestimated subproblem produces a cut, Substituted utility loss (not following WEM) and the server resolves the master problem until the stopping subject to v∈V (35) criterion is met. The result comprises two row-stochastic blocks, ∗ ∗ where V denotes the feasible region of v = [v1 , . . . , vM ]. ZN1 ×Y and ZN2 ×Y , each of size 3 × 3. For example, the The MP iteratively estimates the values of v and checks the constraint comparing x1 and x2 at output y1 is feasibility/optimality of v in the subsequent stage, adding cuts ˜ 1 , x2 ) z2,1 , z1,1 ≤ exp ϵpt d(x to V to refine the solution of MP. Stage 2 - Subproblems: Each subproblem subm validates the with the reverse-direction constraint also included when that guessed values of vm generated by the MP, subject to its own pair is enforced. h

Not following WEM

Local reporting. The server returns each user’s mechanism block with its row order. User 1 privately selects the row corresponding to x2 , and User 2 selects the row corresponding to x5 . Each user samples one yk according to the probabilities in that row and sends only the sampled report for the downstream task. Thus, a 3 × 3 block describes possible reports for all three locations in a neighborhood; the row selection itself remains local to the user. V. P ERFORMANCE E VALUATION In this section, we empirically evaluate the effectiveness of our proposed approximation-based distributed LP framework1 .

mean ± standard deviation. Unless otherwise specified, each user has a neighborhood of size |Nm | = 20, and the low-rank SVD surrogate uses rank r = 5. Computation resource. Our tests were conducted on a machine equipped with an Intel Core i9-13900F processor, featuring 24 cores with a base clock speed ranging from 2.00 GHz to a maximum of 5.60 GHz. The system was configured with 32 GB of DDR5 RAM (32 GB, 4800 MHz) and included an NVIDIA GeForce RTX 4090 graphics card with 24 GB of GDDR6X VRAM. We use the MATLAB optimization toolbox “linearprog” to solve the LP problems.

A. Matrix Approximation Evaluation Datasets. The geographic boundaries and dataset statistics for Comparison of Different Coefficient Approximation Metheach city are summarized below: ods. We evaluated four surrogate families in Table II and • Rome, Italy, bounded by (lat = 41.66, lon = 12.24) in the found that low-rank SVD offers the best privacy-accuracy southwest and (lat = 42.10, lon = 12.81) in the northeast, trade-off, consistently yielding low reconstruction error and including 43,160 nodes and 89,739 edges. low violation ratios under small privacy budgets. Based on this • NYC, USA, bounded by (lat = 40.50, lon = -74.25) in the comparison, our framework adopts the low-rank SVD surrogate southwest and (lat = 40.91, lon = -73.70) in the northeast, as the default approximation model because it provides the best including 55,353 nodes and 139,638 edges. balance between accuracy, scalability, and privacy robustness. • London, UK, bounded by (lat = 51.29, lon = -0.51) in the Detailed comparative results are reported in Section V-A. southwest and (lat = 51.69, lon = 0.28) in the northeast, Performance metrics. In this section, we evaluate two key including 12,820 nodes and 299,524 edges. bd 2 , A bd aspects of the approximated matrices, A Nm ×Y , and Nm For each city, we do not rely solely on the full metropolitan b c ANm ×Y : (1) the privacy loss (formally defined in Eq. (1) in graph but also construct a family of smaller, axis-aligned Definition 1), and (2) the approximation accuracy. subregions designed to probe the robustness of our framework To quantify the overall approximation accuracy, we use under varying spatial extents and network densities. In Rome, the relative Frobenius norm error (or RelErr for simplicity) we systematically partition the area spanning longitudes 12.24- between the original matrix A and its approximation Â: 12.81 and latitudes 41.66-42.10 into a set of contiguous r 2 Pn Pm rectangular windows, yielding induced subgraphs with approxA − Â ij ij i=1 j=1 ∥A − Â∥F imately 3,000 to 15,000 nodes. These regions are concentrated qP RelErr = = , (37) Pm n ∥A∥F around the urban core, where the underlying road network is 2 i=1 j=1 Aij dense and highly irregular. For New York City, we select several windows within longitudes -74.25 to -73.70 and where ∥ · ∥F denotes the Frobenius norm. This metric measures latitudes 40.50 to 40.91, covering lower and mid-Manhattan, the overall deviation of the approximation from the original as well as adjacent neighborhoods. This produces subgraphs matrix as a fraction of the original matrix’s magnitude. A lower with approximately 6,000 to 29,000 nodes. Unlike Rome, value indicates a closer match, while a higher value reflects these subgraphs inherit NYC’s strongly anisotropic, grid-like larger deviations. Normalization by ∥A∥F makes the measure structure. For London, we focus on central districts within scale-independent, allowing fair comparisons across matrices longitudes -0.51 to 0.28 and latitudes 51.29 to 51.69, forming with different magnitudes. d For the distance matrices Adx six moderately sized rectangles that encompass Westminster, m and ANm ×Y , approximation the City of London, and neighboring boroughs. Across all accuracy is particularly critical because it directly determines these cities, these carefully chosen subregions provide compact the mDP indistinguishability constraints. If an approximated ˆ i , xj ) exceeds the true distance dx ,x , the yet representative benchmarks, allowing us to evaluate the distance value d(x i j ˆ i ,xj ) effect of our obfuscation and reconstruction methods across term eϵd(x becomes larger, thereby relaxing the constraint heterogeneous spatial environments. zi,k − eϵdxi ,xj zj,k ≤ 0, (38) Experimental and privacy-auditing settings. We evaluate overall privacy budgets ϵ ∈ {4, 7, 10} km−1 . For each coef- which can lead to solutions that violate the intended privacy ficient type t ∈ {dx, dy, cy}, the upload allocation ϵt is set guarantee under the true distances. To capture this effect, we as from the maximum privacy loss observed in its audit; the define the violation ratio n o ˆ i , xj ) > dx ,x remaining budget is ϵpt = ϵdx + ϵdy + ϵcy . We set the audit (i, j)|d(x i j threshold to δ = 10−4 and use Sℓ = 100 audit samples. Each ViolRatio = , (39) n×n experiment is repeated 10 times, and results are reported as where | · | denotes set cardinality and n × n is the total number 1 Artifact is available at https://anonymous.4open.science/r/DistA-LP-D2F0 of distance entries. A lower violation ratio reflects a more

TABLE II. Performance of coefficient approximation across the three datasets (mean ± 1.96×standard error). Distance matrix Ãdx m RelErr ViolRatio

Distance matrix ÃdNm ×Y RelErr ViolRatio

Utility loss matrix ÃcNm ×Y PL RelErr

Methods

2D Gaussian 2D Polynomial RBF Low-Rank SVD

1.34±0.93 5.36±3.27 8.02±7.15 0.95±0.16

0.60±0.07 168.64±75.42 0.57±0.03 0.07±0.00

0.38±0.04 0.39±0.14 0.41±0.03 0.13±0.01

0.50±0.09 257.21±219.4 0.54±0.04 0.04±0.00

0.44±0.06 0.45±0.16 0.41±0.04 0.22±0.01

1.93±1.91 6.42±4.64 8.06±6.87 1.25±0.23

0.60±0.05 53.83±49.01 0.70±0.02 0.03±0.00

2D Gaussian 2D Polynomial RBF Low-Rank SVD

0.96±0.28 3.24±2.27 3.86±0.90 0.91±0.09

0.85±0.66 220.31±123.4 0.52±0.01 0.07±0.00

London road map 0.42±0.07 0.67±0.29 0.48±0.03 0.43±0.13 2.32±1.15 189.81±116.53 0.43±0.03 3.52±1.14 0.54±0.02 0.13±0.01 0.64±0.12 0.04±0.00

0.45±0.06 0.45±0.17 0.38±0.03 0.20±0.02

1.00±0.56 3.04±0.71 3.36±1.10 0.92±0.58

0.48±0.02 161.41±81.89 0.63±0.01 0.03±0.00

2D Gaussian 2D Polynomial RBF Low-Rank SVD

0.52±0.14 2.43±0.49 3.13±2.12 1.02±0.20

0.66±0.03 253.83±141.4 0.64±0.03 0.08±0.01

New York City road map 0.39±0.02 0.51±0.23 0.66±0.07 0.42±0.22 2.07±0.29 219.51±61.59 0.45±0.03 4.19±1.91 0.64±0.03 0.09±0.01 0.63±0.10 0.06±0.01

0.44±0.08 0.50±0.18 0.42±0.08 0.15±0.02

0.54±0.28 3.32±2.63 3.94±2.87 1.03±0.36

0.54±0.00 48.78±31.61 0.67±0.02 0.03±0.00

Rome road map 1.11±0.66 5.59±3.04 9.26±8.62 0.76±0.13

conservative approximation that preserves or strengthens the • Pre-defined noise distribution, including Exponenoriginal privacy constraints. tial Mechanism (EM) [21], which samples y with Pr[M(x) = y] ∝ exp − ϵd2 (x, y) , where d2 is the Experimental results. Table II compares RelErr, violation ℓ2 distance, and Planar Laplace (Laplace) [6], which ratio, and privacy loss of different approximation methods achieves geo-indistinguishability by adding 2D Laplace across the three datasets. The performance patterns observed in 2 noise, Pr[M(x) = y] = 2πϵ m exp − ϵd2 (x, y) . On Table II are consistent with the modeling/complexity trade-offs discrete grids, Laplace uses an effective ϵ′ ≤ ϵ to maintain summarized in Table VII. Particularly, Low-rank SVD captures mDP. the global structure of the distance/utility matrices with only • Optimization-based methods, including LP [13] and k(m + n) parameters and inherits an implicit denoising effect PAnDA [4]. LP optimizes a mechanism that minimizes from truncation, so its parameter perturbation concentrates on expected loss subject to mDP constraints, which scales a small, well-conditioned set, yielding both low reconstruction poorly over X . PAnDA uses a two-phase framework in error and low violation ratios in Table II. 2D Gaussian offers which each user selects a small anchor set, allowing the a compact, smooth parametric surface; when the underlying server to solve a compact LP on a reduced domain. We matrices are approximately smooth and unimodal within local partition the target region into a coarse 8 × 8 grid and neighborhoods, its few parameters perturb stably, producing apply the LP method (low computation time but inaccurate moderate errors and moderate violation ratios (Table II), making utility loss), denoted as LP. it a reasonable backup where SVD is impractical. RBF has high • Hybrid Methods, including COPT [2], which runs LP expressivity but a much larger parameter set; privacy noise must on a small, utility-sensitive subset of outputs and uses be spread over many coefficients (and often asymmetrically EM elsewhere—cutting cost while keeping mDP, and to curb violations), which amplifies effective perturbation and Bayesian Remapping (RMP) [26], a post-process of an explains the higher privacy loss and violation ratios observed initial output y to ỹ = arg min ỹ∈Y Ex∼Pr[x|y] [L(x, ỹ)], in Table II despite decent fit on certain regions. Finally, 2D where Pr[x | y] follows from prior px and mechanism Z. Polynomial surrogates can be numerically unstable and poorly RMP improves utility without weakening mDP. capture non-polynomial curvature, especially near boundaries; under perturbation, these instabilities magnify, leading to the Table IV reports the expected utility loss across perturbation very large relative errors and no compensating safety gains seen mechanisms on the Rome, London, and New York City roadin Table II. Overall, Table 1’s complexity/regularity arguments map datasets, with values scaled by 10,000. Across all privacy justify the Table II outcome: SVD provides the best privacybudgets ϵ ∈ {4, 7, 10}km−1 and dataset sizes (N = 2, 000, utility-safety balance, with Gaussian serving as a compact 4, 000, 6, 000), D ISTA-LP consistently achieves the lowest loss and robust runner-up, while RBF and Polynomial methods are with narrow confidence bands, outperforming predefined-noise dominated once privacy noise is accounted for. baselines (EM, Laplace) and optimization/hybrid baselines (RMP, LP, PAnDA). B. Utility Loss Averaged over the three road-network datasets and grid sizes, We evaluate utility using the expected distance between the true location x and the reported location M(x) under the same DistA-LP reduces expected utility loss by 50% compared to EM/Laplace and centralized LP, and achieves about a 49% ℓp -norm. Formally, we compute Ex∼p(x),y∼M(x) [dp (x, y)]. We compare our approach against the following representa- reduction relative to PAnDA. Compared to predefined-noise tive methods based on mDP: mechanisms (EM, Laplace), which calibrate noise only to the

TABLE III. Violation ratio (Mean±1.96×standard error). Number of records Privacy budget (km−1 )

2,000 ϵ = 4.0

4,000

ϵ = 7.0

ϵ = 10.0

ϵ = 4.0

6,000

ϵ = 7.0

ϵ = 10.0

ϵ = 4.0

ϵ = 7.0

ϵ = 10.0

0.00±0.00 0.00±0.00 0.00±0.00 ————————0.57±0.04 0.00±0.00 0.00±0.00

Pre-defined EM Noise Distribution Laplace Hybrid RMP Method COPT Optimization LP Based LP-A Methods PAnDA D ISTA-LP†

0.00±0.00 0.00±0.00 0.00±0.00 ————————0.59±0.04 0.00±0.00 0.00±0.00

Rome road map 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 ————————————————0.59±0.0436 0.57±0.04 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00

Pre-defined EM Noise Distribution Laplace Hybrid RMP Method COPT Optimization LP Based LP-A Methods PAnDA D ISTA-LP†

0.00±0.00 0.00±0.00 0.00±0.00 ————————0.50±0.03 0.00±0.00 0.00±0.00

London road map 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 ————————————————0.50±0.03 0.50±0.03 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00

0.00±0.00 0.00±0.00 0.00±0.00 ————————0.50±0.03 0.00±0.00 0.00±0.00

Pre-defined EM Noise Distribution Laplace Hybrid RMP Method COPT Optimization LP Based LP-A Methods PAnDA † D ISTA-LP

0.00±0.00 0.00±0.00 0.00±0.00 ————————0.56±0.01 0.00±0.00 0.00±0.00

New York City road map 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 ————————————————————————0.56±0.01 0.56±0.01 0.56±0.01 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00 0.00±0.00

0.00±0.00 0.00±0.00 0.00±0.00 ————————0.56±0.01 0.00±0.00 0.00±0.00

perturbation magnitude and treat all directions in the metric space isotropically, DistA-LP explicitly optimizes a utilityloss objective over the discretized domain, so it can allocate probability mass toward utility-favorable outputs while still satisfying mDP constraints. Relative to hybrid methods such as COPT and RMP, which either optimize only on a small subset of outputs or improve utility through a post-processing remap of an existing mechanism, DistA-LP optimizes the full perturbation matrix jointly via the WEM+LP formulation, allowing it to adjust all entries coherently instead of partially correcting a predefined mechanism. Finally, compared to centralized LP optimization, DistA-LP leverages surrogate-based coefficient approximation to obtain conservative yet tight distance/utility matrices that exhibit near-zero violation ratios (as reported in Table III) in the Appendix; this avoids the over-conservative noise that arises when centralized LP accumulates large distance overestimates, thereby yielding lower distortion at similar privacy budgets. C. mDP Violation Ratio Table III reports the empirical mDP violation ratio, defined as the fraction of tested mDP constraints whose privacy loss exceeds the target budget ϵ. A lower violation ratio indicates stronger empirical compliance with the required mDP constraints. Across all three road-network datasets, privacy budgets, and domain sizes, D ISTA-LP achieves 0.00 ± 0.00 violation ratio, matching the predefined mechanisms and postprocessing baselines that preserve mDP by construction. In contrast, LP-A exhibits substantial violation ratios, around 0.50– 0.59, because its coarse approximation can relax distance-based mDP constraints and produce mechanisms that violate the intended privacy bound. These results show that the proposed

surrogate-based approximation and conservative reconstruction steps preserve empirical mDP validity while enabling scalable optimization. D. Computation Efficiency Table V reports the end-to-end computation time for each perturbation mechanism across privacy budgets ϵ ∈ {4, 7, 10}km−1 and dataset sizes (|X | = 2, 000, 4, 000, 6, 000) on the road-network datasets. Predefined-noise baselines (EM, Laplace) are extremely fast (≤ 0.05 time units) and essentially insensitive to ϵ or |X |, since they avoid solving large optimization problems. Hybrid methods (RMP, PAnDA) incur modest overhead (roughly 0.03–0.6 seconds), remaining lightweight compared to full LP. In contrast, centralized LP/COPT often reach the 1,800-second timeout and, when they do finish, exhibit a steep growth in runtime—from about 15.5–18.0 seconds at |X | = 2,000 to ≈ 30.2–37.7 seconds at |X | = 4,000 and ≈ 45.1–46.9 at |X | = 6,000, reflecting their poor scalability. DistA-LP lies between these extremes: it is somewhat more expensive than purely closedform or hybrid baselines but remains in a practical range and shows substantially more stable scaling than centralized LP. Overall, the table illustrates that DistA-LP achieves a favorable trade-off between computation cost and utility, dramatically improving scalability over centralized LP-based mechanisms while retaining their optimization benefits. VI. R ELATED W ORK Since the introduction of mDP by [5], substantial research efforts have focused on developing data perturbation mechanisms. These mechanisms generally fall into two categories:

TABLE IV. Utility loss across different perturbation methods (10,000 meters) Mean±1.96×standard error; the algorithm without getting the results, we label its utility by “——–”. Number of records −1

Privacy budget (km

2,000 )

ϵ = 4.0

4,000 ϵ = 4.0

6,000

ϵ = 7.0

ϵ = 10.0

ϵ = 7.0

ϵ = 10.0

ϵ = 4.0

ϵ = 7.0

ϵ = 10.0

7.12±2.51 7.30±2.66 6.79±2.54 ————————6.73±2.55 6.43±2.60 2.88±0.21

7.10±2.50 7.30±2.66 6.78±2.54 ————————6.72±2.55 6.43±2.60 2.88±0.21

7.98±3.35 8.01±3.35 7.73±3.37 ————————7.63±3.32 7.39±3.25 3.14±0.26

7.96±3.34 8.01±3.35 7.71±3.36 ————————7.60±3.32 7.39±3.25 3.14±0.26

7.94±3.34 8.01±3.35 7.69±3.35 ————————7.60±3.32 7.39±3.25 3.14±0.26

Pre-defined EM Noise Distribution Laplace Hybrid RMP Method COPT Optimization LP Based LP-A Methods PAnDA † D ISTA-LP

5.63±0.95 5.68±0.95 5.30±0.91 ————————5.19±0.91 4.95±1.03 3.47±0.24

5.60±0.95 5.68±0.95 5.26±0.91 ————————5.15±0.91 4.95±1.03 3.47±0.24

Rome road map 5.57±0.95 7.16±2.53 5.68±0.95 7.30±2.66 5.23±0.92 6.81±2.54 ————————————————5.14±0.91 6.76±2.55 4.95±1.03 6.43±2.60 3.47±0.24 2.88±0.21

Pre-defined EM Noise Distribution Laplace Hybrid RMP Method COPT Optimization LP Based LP-A Methods PAnDA † D ISTA-LP

6.61±2.03 6.74±2.00 6.22±1.99 ————————6.13±2.01 6.03±2.17 2.74±0.20

6.54±2.04 6.74±2.01 6.17±2.00 ————————6.08±2.02 6.02±2.17 2.74±0.21

London road map 6.50±2.05 6.57±1.85 6.74±2.01 6.66±1.86 6.14±2.01 6.30±1.87 ————————————————6.07±2.02 6.22±1.88 6.02±2.17 5.93±1.80 2.74±0.20 3.79±0.08

6.53±1.84 6.66±1.86 6.27±1.87 ————————6.18±1.88 5.93±1.80 3.79±0.08

6.50±1.84 6.66±1.86 6.25±1.87 ————————6.17±1.88 5.93±1.80 3.79±0.08

6.47±1.90 6.51±1.90 6.22±1.92 ————————6.14±1.92 5.88±1.99 3.99±0.11

6.44±1.91 6.51±1.90 6.20±1.92 ————————6.11±1.92 5.88±1.99 3.99±0.11

6.42±1.91 6.51±1.90 6.18±1.92 ————————6.10±1.92 5.88±1.99 3.99±0.11

Pre-defined EM Noise Distribution Laplace Hybrid RMP Method COPT Optimization LP Based LP-A Methods PAnDA D ISTA-LP†

5.64±0.91 5.80±0.92 5.22±0.96 ————————5.12±0.92 5.12±1.20 1.60±0.09

5.58±0.92 5.80±0.92 5.18±0.95 ————————5.05±0.92 5.12±1.20 1.58±0.06

New York City road map 5.54±0.93 5.66±0.85 5.62±0.85 5.80±0.92 5.76±0.87 5.77±0.87 5.15±0.95 5.36±0.84 5.34±0.84 ————————————————————————5.04±0.92 5.25±0.83 5.22±0.84 5.12±1.21 4.78±0.73 4.78±0.73 1.54±0.03 2.25±0.16 2.25±0.17

5.59±0.85 5.77±0.87 5.32±0.84 ————————5.21±0.84 4.78±0.73 2.25±0.15

4.98±0.81 5.14±0.87 4.74±0.75 ————————4.66±0.74 4.89±0.98 2.47±0.09

4.93±0.79 5.14±0.87 4.72±0.75 ————————4.62±0.73 4.89±0.98 2.46±0.11

4.91±0.79 5.14±0.87 4.70±0.74 ————————4.61±0.73 4.89±0.98 2.45±0.09

predefined noise distribution mechanisms and optimizationbased mechanisms.

many recent efforts have prioritized optimization-based mechanisms, which minimize data utility loss by explicitly modeling both the magnitude and directional impact of each candidate perturbation.

Predefined Noise Distribution Mechanisms. The Laplace mechanism is a classic perturbation mechanism, originally proposed for DP [27] and later adapted to mDP by Optimization (or LP)-based Mechanisms. Since the seminal Chatzikokolakis et al. [5]. Instead of scaling noise to global work by Bordenabe et al. [13], which introduced the use of LP sensitivity, it uses the metric distance dX between records to optimize perturbation mechanisms under mDP, significant and sets the noise scale to dϵX , so that outputs for nearby effort has been devoted to developing LP-based mechanisms. records remain statistically indistinguishable. This adaptation These approaches generally aim to minimize utility loss by underlies Geo-Indistinguishability [6] for location privacy: explicitly considering both the magnitude and direction of early work focused on sporadic locations [5], [6], while the loss for each possible perturbation. Such mechanisms later extensions considered continuous tracking, real-time are typically formalized as stochastic matrices, with the LP trajectories, and personalized budgets [28], [29], [30], [31]. objective minimizing the expected utility loss subject to Despite its efficiency [3], additive Laplace noise is poorly mDP constraints. However, these formulations involve O(N 2 ) suited to structured domains (e.g., graphs, embeddings), decision variables, making them theoretically optimal but where it can distort structure or ignore density. For word computationally prohibitive for large domains. Many existing embeddings, perturbation may produce invalid vectors, implementations (e.g., [13], [44], [7], [22]) are limited to requiring nearest-neighbor projection [32], [33], which in turn relatively small domains with N ≤ 100. Although Bordenabe assumes the embedding space itself is non-sensitive and can et al. [13] proposed using spanners to approximate dX via introduce privacy risks during neighbor retrieval. shortest path distances in a graph, the scalability remains EM [3], [2] addresses these limitations by selecting outputs limited, for example, solving the LP takes over 1,800 seconds from discrete sets via utility-driven quality scores, bypassing when N = 400 [2]. continuous perturbation. This proves particularly effective Recent research has prioritized improving the scalability in location privacy (trajectory protection [21], [7]; spatial of LP-based methods. Ahuja et al. [45] proposed a hierarcrowdsourcing [34], [35], [36], [37]), textual data (word chical index structure to protect location data under Geosubstitution [38], [39]; embeddings [3], [40]), and graph data Indistinguishability (Geo-Ind), while others leveraged decompo(edge preservation [41]; query responses [42], [43]). While EM sition techniques, such as Dantzig-Wolfe [12], [46] and Benders offers greater flexibility and often superior utility compared to decomposition [14], to partition large-scale LP problems into the Laplace mechanism, it bases selection on perturbation tractable subproblems. Hybrid frameworks combining LP with magnitude dx,y , which may not fully capture directional probabilistic mechanisms have also emerged, adapting noise variations in utility loss across the output space. Therefore, using local utility/sensitivity. For example, Imola et al. [2]

TABLE V. Computation time (Seconds) Mean±1.96×standard error. Number of records Privacy budget (km−1 )

2,000 ϵ = 4.0

4,000

ϵ = 7.0

ϵ = 10.0

ϵ = 4.0

6,000

ϵ = 7.0

ϵ = 10.0

ϵ = 4.0

ϵ = 7.0

ϵ = 10.0

≤ 0.05 ≤ 0.05 0.26±0.01 ≥ 1, 800 ≥ 1, 800 23.03±0.67 0.163±0.070 4.55±0.11

≤ 0.05 ≤ 0.05 0.27±0.01 ≥ 1, 800 ≥ 1, 800 28.67±0.95 0.161±0.065 4.54±0.19

≤ 0.05 ≤ 0.05 0.45±0.00 ≥ 1, 800 ≥ 1, 800 32.24±0.41 0.33±0.30 2.75±0.21

≤ 0.05 ≤ 0.05 0.45±0.00 ≥ 1, 800 ≥ 1, 800 33.33±0.53 0.29±0.22 2.74±0.36

≤ 0.05 ≤ 0.05 0.58±0.04 ≥ 1, 800 ≥ 1, 800 36.95±1.21 0.29±0.21 3.71±0.52

Pre-defined EM Noise Distribution Laplace Hybrid RMP Method COPT Optimization LP Based LP-A Methods PAnDA D ISTA-LP†

≤ 0.05 ≤ 0.05 ≤ 0.05 ≥ 1, 800 ≥ 1, 800 11.27±0.12 0.083±0.014 2.03±0.16

≤ 0.05 ≤ 0.05 ≤ 0.05 ≥ 1, 800 ≥ 1, 800 11.79±0.86 0.083±0.016 2.07±0.20

Rome road map ≤ 0.05 ≤ 0.05 ≤ 0.05 ≤ 0.05 ≤ 0.05 0.25±0.00 ≥ 1, 800 ≥ 1, 800 ≥ 1, 800 ≥ 1, 800 13.69±0.54 21.99±0.29 0.082±0.015 0.157±0.057 2.15±0.33 4.68±0.36

Pre-defined EM Noise Distribution Laplace Hybrid RMP Method COPT Optimization LP Based LP-A Methods PAnDA D ISTA-LP†

≤ 0.05 ≤ 0.05 0.03±0.00 ≥ 1, 800 ≥ 1, 800 41.8±0.51 0.09±0.01 2.08±0.05

≤ 0.05 ≤ 0.05 0.03±0.00 ≥ 1, 800 ≥ 1, 800 42.9±1.99 0.09±0.00 3.89±2.76

London road map ≤ 0.05 ≤ 0.05 ≤ 0.05 ≤ 0.05 0.03±0.00 0.20±0.01 ≥ 1, 800 ≥ 1, 800 ≥ 1, 800 ≥ 1, 800 42.1±0.52 83.9±2.06 0.09±0.01 0.14±0.02 3.21±0.12 0.51±0.06

≤ 0.05 ≤ 0.05 0.18±0.00 ≥ 1, 800 ≥ 1, 800 80.7±0.94 0.14±0.02 0.77±0.15

≤ 0.05 ≤ 0.05 0.22±0.01 ≥ 1, 800 ≥ 1, 800 87.2±2.87 0.15±0.03 0.67±0.01

≤ 0.05 ≤ 0.05 0.47±0.04 ≥ 1, 800 ≥ 1, 800 127.6±9.8 0.24±0.04 2.22±0.23

≤ 0.05 ≤ 0.05 0.58±0.07 ≥ 1, 800 ≥ 1, 800 145.8±13.7 0.24±0.06 2.31±0.37

≤ 0.05 ≤ 0.05 0.55±0.03 ≥ 1, 800 ≥ 1, 800 144.1±2.6 0.24±0.03 2.23±0.28

Pre-defined EM Noise Distribution Laplace Hybrid RMP Method COPT Optimization LP Based LP-A Methods PAnDA † D ISTA-LP

≤ 0.05 ≤ 0.05 ≤ 0.05 ≥ 1, 800 ≥ 1, 800 15.52±0.50 0.089±0.015 0.92±0.02

≤ 0.05 ≤ 0.05 ≤ 0.05 ≥ 1, 800 ≥ 1, 800 15.84±0.52 0.090±0.022 0.94±0.05

New York City road map ≤ 0.05 ≤ 0.05 ≤ 0.05 ≤ 0.05 ≤ 0.05 ≤ 0.05 ≤ 0.05 0.18±0.00 0.23±0.01 ≥ 1, 800 ≥ 1, 800 ≥ 1, 800 ≥ 1, 800 ≥ 1, 800 ≥ 1, 800 18.03±1.13 30.19±0.52 37.70±1.82 0.089±0.018 0.13±0.017 0.14±0.0063 0.93±0.05 18.28±0.07 18.57±0.18

≤ 0.05 ≤ 0.05 0.22±0.01 ≥ 1, 800 ≥ 1, 800 33.62±0.39 0.15±0.021 18.70±0.18

≤ 0.05 ≤ 0.05 0.42±0.01 ≥ 1, 800 ≥ 1, 800 45.09±1.23 0.24±0.078 2.77±0.31

≤ 0.05 ≤ 0.05 0.45±0.00 ≥ 1, 800 ≥ 1, 800 46.00±1.05 0.23±0.073 2.81±0.37

≤ 0.05 ≤ 0.05 0.53±0.01 ≥ 1, 800 ≥ 1, 800 46.92±1.22 0.23±0.084 2.80±0.45

integrated the weighted EM with selective LP optimization on sensitive inputs. Despite these advances, these state-of-theart methods still only handle domains with N ≤ 1,000. A more recent work by Qiu et al. [15] introduces a localityaware optimization approach, restricting computation to spatial neighborhoods. While scaling to N = 1,600, it does so at the cost of strict compliance with mDP, as the use of neighborhoodspecific parameters introduces potential information leakage. VII. C ONCLUSIONS We introduced D ISTA-LP, a privacy-preserving distributed LP framework for scalable mDP optimization with explicit upload-stage privacy accounting. D ISTA-LP replaces raw user-specific LP coefficients with noisy surrogate parameters, composes upload leakage with the downstream perturbation mechanism for end-to-end (ϵ, δ) guarantees, and uses privacyfree post-processing plus Benders-style decomposition for scalable, geometrically consistent optimization. Experiments on Rome, London, and NYC road-network datasets show that D ISTA-LP achieves lower utility loss, near-zero reported mDP violations, and substantially faster runtime than monolithic LP and hybrid baselines. These results demonstrate that surrogate parameterization and probabilistic auditing can protect uploadstage privacy without sacrificing utility or formal guarantees. Promising avenues include deriving tighter analytic bounds for asymmetric parameter noise to complement the PmDP auditor; exploring richer surrogate families (e.g., kernelized or learned surrogates) with privacy-aware regularization; extending the framework beyond geospatial domains to text, graphs, and multimodal metrics; handling adversarial or colluding participants in stronger threat models; and integrating adaptive privacy budgeting across time or tasks. Finally, we believe

D ISTA-LP provides a principled foundation for scalable, auditable mDP optimization in distributed settings. E THICAL C ONSIDERATIONS We reviewed the venue’s ethics guidance, submission instructions, and ethics-document requirements. This work was conducted as a simulation-based study using public infrastructurelevel data, and we believe the research process poses minimal direct risk to individuals. At the same time, because the paper concerns privacy mechanisms for sensitive records, especially location data, we explicitly consider potential downstream risks and limitations. a) Stakeholders: The primary stakeholders include: (i) researchers studying metric differential privacy (mDP), distributed optimization, and privacy auditing; (ii) practitioners and platform operators who may deploy mDP mechanisms in navigation, mobility analytics, spatial crowdsourcing, or related services; (iii) end users and data subjects whose locations or other sensitive records may be protected by such mechanisms; (iv) data and infrastructure providers such as OpenStreetMap contributors; and (v) society at large, including communities that may face heightened risks from location tracking, surveillance, or privacy-invasive analytics. The research process itself has limited direct impact on these stakeholders, while publication and deployment may influence how distributed mDP mechanisms are designed, audited, and used in practice. b) Impacts of the research process and publication: Regarding the research process, our evaluation uses publicly available OpenStreetMap road-network graphs and synthetic perturbation experiments. We do not collect, process, or attempt to infer personal or identifiable user-level trajectories, and we

do not interact with live users or deployed systems. Thus, and giving practitioners concrete guidance for designing the risk of direct tangible harms or rights-based harms to and auditing privacy-preserving distributed LP mechanisms. individuals during the research process is minimal. We believe that making the assumptions, limitations, and Regarding the publication and potential deployment of the safeguards explicit is preferable to leaving upload-stage leakage results, the intended impact is positive. DISTA-LP is designed unexamined in future distributed privacy systems. to make LP-based mDP optimization more scalable while R EFERENCES explicitly accounting for privacy leakage from distributed coefficient uploads. In particular, the framework replaces raw [1] J. C. Duchi, M. I. Jordan, and M. J. Wainwright, “Local privacy and user-specific LP coefficient blocks with perturbed surrogate statistical minimax rates,” in 2013 IEEE 54th Annual Symposium on Foundations of Computer Science, 2013, pp. 429–438. parameters, uses probabilistic mDP auditing to quantify upload[2] J. Imola, S. Kasiviswanathan, S. White, A. Aggarwal, and N. Teissier, stage leakage, and composes this leakage with the downstream “Balancing utility and scalability in metric differential privacy,” in Proc. perturbation mechanism under an end-to-end privacy budget. of UAI 2022, 2022. [3] R. S. Carvalho, T. Vasiloudis, and O. Feyisetan, “Tem: High utility These contributions can benefit researchers by clarifying a metric differential privacy on text,” ArXiv, vol. abs/2107.07928, 2021. previously under-accounted leakage channel, practitioners by [Online]. Available: https://api.semanticscholar.org/CorpusID:236034456 providing an auditable design pattern for distributed mDP [4] R. Liu and C. Qiu, “Panda: Rethinking metric differential privacy optimization at scale with anchor-based approximation,” in Proc. of optimization, and data subjects by strengthening protections ACM CCS, 2025. in settings where local distance or utility structures may [5] K. Chatzikokolakis, M. E. Andrés, N. E. Bordenabe, and C. Palamidessi, themselves reveal sensitive information. “Broadening the scope of differential privacy using metrics,” in Proc. of Privacy Enhancing Technologies, 2013, pp. 82–102. At the same time, several risks remain. First, the framework [6] M. E. Andrés, N. E. Bordenabe, K. Chatzikokolakis, and C. Palamidessi, could be misused to justify weak privacy configurations, such “Geo-indistinguishability: Differential privacy for location-based systems,” as overly large ϵ or δ values, poorly chosen metrics, or in Proc. of ACM CCS, 2013, pp. 901–914. [7] L. Yu, L. Liu, and C. Pu, “Dynamic differential location privacy with insufficient sampling in the PmDP auditor. Second, practitioners personalized error bounds,” in Proc. of IEEE NDSS, 2017. may over-interpret the empirical guarantees outside the stated [8] R. Shokri, “Privacy games: Optimal user-centric data obfuscation,” Proc. assumptions, such as applying the method under stronger on Privacy Enhancing Technologies, vol. 2015, no. 2, pp. 299 – 315, 2015. adversaries, colluding users, or domains whose metric does [9] F. Ahmed, A. X. Liu, and R. Jin, “Social graph publishing with privacy not reflect meaningful privacy harms. Third, because the paper guarantees,” in Proc. of IEEE ICDCS, 2016. highlights leakage from coefficient uploads, careless readers [10] Y. Han et al., “Voice-indistinguishability: Protecting voiceprint in privacypreserving speech data release,” in IEEE ICME, 2020. might focus on this vulnerability without adopting the proposed mitigation. These risks are particularly important for high-risk [11] J.-W. Chen et al., “Perceptual indistinguishability-net: Facial image obfuscation with manipulable semantics,” in Proc. of CVPR, 2021. location-based applications and for populations that may be [12] C. Qiu, A. C. Squicciarini, C. Pang, N. Wang, and B. Wu, “Location disproportionately affected by surveillance or privacy failures. privacy protection in vehicle-based spatial crowdsourcing via geoindistinguishability,” IEEE Transactions on Mobile Computing, pp. 1–1, c) Mitigations and residual risks: To mitigate these risks, 2022. we clearly document the threat model, including the honest- [13] N. E. Bordenabe, K. Chatzikokolakis, and C. Palamidessi, “Optimal but-curious server assumption, honest users, and the exclusion geo-indistinguishable mechanisms for location privacy,” in Proc. of ACM CCS, 2014. of collusion from the main analysis. We also distinguish exact [14] C. Qiu, “Enhancing scalability of metric differential privacy via secret mDP from probabilistic mDP and emphasize that the uploaddataset partitioning and benders decomposition,” in Proc. of IJCAI, 2024. stage privacy budget must be explicitly allocated and composed [15] C. Qiu, R. Liu, P. Pappachan, A. Squicciarini, and X. Xie, “Time-efficient locally relevant geo-location privacy protection,” in Privacy Enhancing with the downstream perturbation budget. We recommend Technologies Symposium (PETS), 2025. that practitioners use conservative privacy budgets, validate [16] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed. Baltimore, metric choices with respect to the relevant application and MD: Johns Hopkins University Press, 2013. affected population, perform adequate privacy auditing before [17] G. A. F. Seber and A. J. Lee, Linear Regression Analysis, 2nd ed. NJ: Wiley, 2012. deployment, and avoid deploying the method in high-stakes [18] Hoboken, C. M. Bishop, Pattern Recognition and Machine Learning. New York, settings without additional legal, ethical, and domain-specific NY: Springer, 2006. review. Our released artifacts are limited to code and public [19] M. D. Buhmann, Radial Basis Functions: Theory and Implementations. UK: Cambridge University Press, 2003. or derived experimental data needed for reproducibility; they [20] Cambridge, “openstreetmap,” https://www.openstreetmap.org/, 2020, accessed: 2020do not include personal trajectories or sensitive raw user data. 04-07. Despite these mitigations, residual risk remains that third parties [21] K. Chatzikokolakis, C. Palamidessi, and M. Stronati, “Constructing elastic distinguishability metrics for location privacy,” Privacy Enhancing may ignore the assumptions, select weak parameters, or deploy Technologies (PoPETs), vol. 2015, p. 156–170, 2015. the method in inappropriate contexts. [22] L. Wang, D. Yang, X. Han, T. Wang, D. Zhang, and X. Ma, “Location privacy-preserving task allocation for mobile crowdsensing with d) Decision to conduct and publish the research: We differential geo-obfuscation,” in Proc. of ACM WWW, 2017, pp. 627–636. judged that conducting this research is ethically justified [23] H. To, G. Ghinita, L. Fan, and C. Shahabi, “Differentially private because it uses public infrastructure-level data and synthetic location protection for worker datasets in spatial crowdsourcing,” IEEE Transactions on Mobile Computing, pp. 934–949, 2017. experiments, while addressing an important privacy risk in distributed mDP optimization: leakage from user-specific [24] K. Chatzikokolakis et al., “Constructing elastic distinguishability metrics for location privacy,” Proc. on PETS, 2015. coefficient uploads. Publishing the work provides a net benefit [25] J. A. Nelder and R. Mead, “A simplex method for function minimization,” by enabling community scrutiny, improving reproducibility, The Computer Journal, vol. 7, no. 4, pp. 308–313, 1965.

[26] K. Chatzikokolakis, E. Elsalamouny, and C. Palamidessi, “Efficient utility improvement for location privacy,” Proceedings on Privacy Enhancing Technologies, vol. 2017, 10 2017. [27] C. Dwork, F. McSherry, K. Nissim, and A. Smith, “Calibrating noise to sensitivity in private data analysis,” in Theory of Cryptography. Springer Berlin Heidelberg, 2006, pp. 265–284. [28] K. Chatzikokolakis et al., “A predictive differentially-private mechanism for mobility traces,” in Proc. on PETS, 2014. [29] J. Hua, W. Tong, F. Xu, and S. Zhong, “A geo-indistinguishable location perturbation mechanism for location-based services supporting frequent queries,” IEEE Transactions on Information Forensics and Security, vol. 13, no. 5, pp. 1155–1168, 2018. [30] L. Yu, S. Zhang, Y. Meng, S. Du, Y. Chen, Y. Ren, and H. Zhu, “Privacy-preserving location-based advertising via longitudinal geoindistinguishability,” IEEE Transactions on Mobile Computing, vol. 23, no. 8, pp. 8256–8273, 2024. [31] M. Min, H. Zhu, J. Ding, S. Li, L. Xiao, M. Pan, and Z. Han, “Personalized 3d location privacy protection with differential and distortion geo-perturbation,” IEEE TDSC, vol. 21, no. 4, pp. 3629–3643, 2024. [32] N. Fernandes, M. Dras, and A. McIver, “Author obfuscation using generalised differential privacy,” arXiv preprint arXiv:1805.08866, 2018. [33] O. Feyisetan, T. Diethe, and T. Drake, “Leveraging hierarchical representations for preserving privacy and utility in text,” in 2019 IEEE International Conference on Data Mining (ICDM), nov 2019, pp. 210– 219. [34] M. E. Gursoy, A. Tamersoy, S. Truex, W. Wei, and L. Liu, “ Secure and Utility-Aware Data Collection with Condensed Local Differential Privacy ,” IEEE Transactions on Dependable and Secure Computing, vol. 18, no. 05, pp. 2365–2378, Sep. 2021. [35] B. Niu, Y. Chen, Z. Wang, F. Li, B. Wang, and H. Li, “Eclipse: Preserving differential location privacy against long-term observation attacks,” IEEE Transactions on Mobile Computing, vol. 21, no. 1, pp. 125–138, 2022. [36] Y. Ren, X. Li, Y. Miao, R. H. Deng, J. Weng, S. Ma, and J. Ma, “Distpreserv: Maintaining user distribution for privacy-preserving locationbased services,” IEEE Transactions on Mobile Computing, vol. 22, no. 6, pp. 3287–3302, 2023. [37] X. Dong, T. Zhang, D. Lu, G. Li, Y. Shen, and J. Ma, “Preserving geoindistinguishability of the primary user in dynamic spectrum sharing,” IEEE Transactions on Vehicular Technology, vol. 68, no. 9, pp. 8881– 8892, 2019. [38] S. Wang, Y. Nie, P. Wang, H. Xu, W. Yang, and L. Huang, “Local private ordinal data distribution estimation,” in IEEE INFOCOM 2017-IEEE Conference on Computer Communications. IEEE, 2017, pp. 1–9. [39] O. Feyisetan, B. Balle, T. Drake, and T. Diethe, “Privacy- and utilitypreserving textual analysis via calibrated multivariate perturbations,” in Proc. of ACM WSDM, 2020, p. 178–186. [40] S. Meisenbacher, M. Chevli, and F. Matthes, “1-diffractor: Efficient and utility-preserving text obfuscation leveraging word-level metric differential privacy,” arXiv preprint arXiv:2405.01678, 2024. [41] S. Raskhodnikova and A. Smith, “Lipschitz extensions for node-private graph statistics and the generalized exponential mechanism,” in IEEE Annual Symposium on Foundations of Computer Science (FOCS). IEEE, 2016, pp. 495–504. [42] F. Fioretto, T. W. Mak, and P. Van Hentenryck, “Privacy-preserving obfuscation of critical infrastructure networks,” arXiv preprint arXiv:1905.09778, 2019. [43] P. Kamalaruban, V. Perrier, H. J. Asghar, and M. A. Kaafar, “Not all attributes are created equal: dx-private mechanisms for linear queries,” Proceedings on Privacy Enhancing Technologies, 2020. [44] L. Wang, D. Zhang, D. Yang, B. Lim, and X. Ma, “Differential location privacy for sparse mobile crowdsensing,” 12 2016, pp. 1257–1262. [45] R. Ahuja, G. Ghinita, and C. Shahabi, “A utility-preserving and scalable technique for protecting location data with geo-indistinguishability,” in EDBT, 2019, pp. 217–228. [46] C. Qiu, A. C. Squicciarini, Z. Li, C. Pang, and L. Yan, “Time-efficient geo-obfuscation to protect worker location privacy over road networks in spatial crowdsourcing,” in Proc. of ACM CIKM, 2020.

A PPENDIX This appendix provides supplementary details supporting the main paper, including extended derivations, detailed algorithm design, and complete results omitted from the main text due

to space constraints. In particular, it presents the mathematical details of the proposed DistA-LP framework, privacy-loss accounting procedures, surrogate approximation methods, and additional empirical analyses used to validate scalability, privacy preservation, and utility performance. A. Mathematical Notations TABLE VI. Notation summary Sets & Indices X Y M T

Secret/input space; records xi ∈ X . Output/reported space; yk ∈ Y. Perturbation mechanism; M : X → Y Index set of coefficient-block types; T = {dx, dy, cy} Metric on X ; di,j = dxi ,xj . Number of users; m ∈ {1, . . . , M }. User-m local neighborhood (relevance subset). “Far” set for m handled via WEM structure. Record indices in X . Output index in Y.

d(·, ·) M Nm Hm i, j k

Matrices, Variables, and Parameters |X |×|Y|

Z ∈ R≥0 zi,k Adx m

Ady m Acy m A fΘ e Θ, Θ e = f (·; Θ) e A r U, Σ, V P, Q w vm

P Mechanism matrix; row i is zi⊤ and k zi,k = 1. Pr[yk | xi ]. LP coefficients for indistinguishability over (i, j) ∈ 2 Nm . Distance-based LP coefficients for (i, k) ∈ Nm ×Y. Cost-based LP coefficients for (i, k) ∈ Nm × Y. Generic placeholder for a coefficient matrix above. Coefficient surrogate with parameters Θ. Unperturbed / perturbed surrogate parameters. Reconstructed approximate coefficients on the server. Surrogate rank (SVD/low-rank). Truncated SVD factors: A ≈ UΣV⊤ . Low-rank factors: A ≈ PQ⊤ (uploads may use e Q). e P, WEM weights (master problem). User-m subproblem value in Benders cuts.

Objectives, Constraints, and Privacy Metrics FLP (·) Ω(·) ϵ δ Lxi ,xj (y) ViolRatio RelErr

Local LP over (Nm , Y) with (D, C) and (Ad , Ac ). Feasible region (local/master LP). mDP privacy parameter. Failure probability for probabilistic mDP (PmDP). Pr[M(xi )=y] ln Pr[M(x dxi ,xj (privacy loss). j )=y] Fraction of violations of ϵ-mDP. Relative error vs.true coefficients/utilities.

B. Omitted Proofs Proof of Proposition 1 Proposition 1. (Sequential composition of PmDP, adapted from [4]) Let {Mpt }t∈T be randomized mechanisms with conditionally independent outputs P given the input. If each Mpt satisfies (ϵt , δt )-PmDP and t∈T δtP < 1, thenPthe joint mechanism M = (Mpt )t∈T satisfies ( t∈T ϵt , t∈T δt )PmDP.

Proof. Let z = (zt )t∈T denote the joint output of M. Since the mechanisms are conditionally independent given the input, the privacy loss of the joint mechanism satisfies Pr[M(xi ) = z] PLM xi ,xj (z) = log Pr[M(xj ) = z] X = PL(t) xi ,xj (zt ),

for t ∈ {dx, dy, cy}. For i.i.d.Laplace noise, the privacy loss between two secrets xi , xj satisfies

t∈T p where PL(t) xi ,xj (zt ) denotes the privacy loss of Mt .

By the union bound, for any input pair (xi , xj ), we have " # X M Pr PLxi ,xj (z) > ϵt t∈T (41) i X h (z ) > ϵ . ≤ Pr PL(t) t t xi ,xj

PLxi ,xj (y) = log

∆t ≜ max i̸=j

(42)

Proof of Proposition 2 (Concentration Guarantees) Proposition 2 (Concentration Guarantee for Probabilistic mDP Sampling). If the empirical violation hrate satisfies i p̂Sℓ = ξδ for some constant ξ < 1, then Pr pXℓ2 ≤ δ ≥ 1 − 2 exp(−2Sℓ (1 − ξ)2 δ 2 ). Proof. Suppose the empirical violation rate satisfies p̂S = ξδ for some constant ξ < 1. Our goal is to show that, with high probability, the true violation probability pX 2 is at most δ. First, by Hoeffding’s inequality, for any t > 0: (43)

Let t = (1 − ξ)δ. Then: 2 2

(44)

This implies: 2 2

Pr [|p̂S − pX 2 | ≤ (1 − ξ)δ] ≥ 1 − 2e−2Sℓ (1−ξ) δ .

(45)

Now, under the event that |p̂S − pX 2 | ≤ (1 − ξ)δ, we have: pX 2 ≤ p̂S + (1 − ξ)δ = ξδ + (1 − ξ)δ = δ.

Thus:

2 2

Pr [pX 2 ≤ δ] ≥ 1 − 2e−2Sℓ (1−ξ) δ ,

vec(Atm (xi )) − vec(Atm (xj )) 1 . d(xi , xj )

reported in Fig. 3(a)–(c).

where the last inequality follows from the (ϵt , δt )-PmDP guarantee of each Mpt . P P Therefore, by Definition 3, M satisfies ( t∈T ϵt , t∈T δt )PmDP.

Pr [|p̂S − pX 2 | > (1 − ξ)δ] ≤ 2e−2Sℓ (1−ξ) δ .

vec(Atm (xi )) − vec(Atm (xj )) 1 . bt

In our experiments, we fix ϵ = 4 and apply entrywise perturba t e then compute the induced privacy A, tion Lap 0, ∆ to obtain ϵ e . The resulting privacy-loss distributions are loss PLx ,x A

t∈T

p(y | xi ) ≤ p(y | xj )

To align with ϵ-mDP, i.e., PLxi ,xj (y) ≤ ϵ d(xi , xj ) for all i, j, y, it suffices to set bt = ∆t /ϵ, where the metric (Lipschitz) sensitivity is

t∈T

Pr [|p̂S − pX 2 | > t] ≤ 2e−2Sℓ t .

Laplace Perturbation Baseline for Coefficient Blocks We instantiate a naive additive-noise baseline by releasing Mt (x) = vec Atm (x) + η, η ∼ i.i.d. Lap(0, bt ),

(40)

Averaging over all input pairs yields " # X X 1 M Pr PLxi ,xj (z) > ϵt |X |2 xi ,xj ∈X t∈T h i X X 1 (t) (z ) > ϵ Pr PL ≤ t t xi ,xj |X |2 xi ,xj ∈X t∈T X ≤ δt ,

C. Detailed Algorithm Design

(46) (47)

which completes the proof.

Coefficient Matrix Approximation Joint Optimization of Permutation {Em , Θtm }M m=1 : (1) Parameter fitting. Holding Em fixed, for each surrogate family we fit Θtm by minimizing the weighted entrywise ℓ1 reconstruction loss in (23) using the Nelder–Mead simplex algorithm [25], a robust, derivative-free optimizer wellsuited for nonconvex objectives and nonsmooth losses. This choice avoids instabilities we observed with gradientbased methods under discretization and privacy-induced perturbations. For the Gaussian family, we impose implicit box constraints to ensure physically meaningful spreads (sx , sy ) and nonnegative amplitude (A ≥ 0), avoiding degenerate fits. For the polynomial family, we add mild ridge regularization to improve conditioning and reduce overfitting, i.e., we minimize the ℓ1 reconstruction loss plus ρ∥θ∥22 over polynomial coefficients. For the RBF family, we restrict the bandwidth σ to a bounded range for numerical stability and softly regularize the basis weights αk to prevent overconcentration and preserve smoothness across neighboring centers. (2) Permutation search. With {Θtm } fixed, we refine the local permutation Em via greedy pairwise index swaps (i, j) that decrease the global objective in (23). We accept only swaps that reduce the total weighted ℓ1 loss, ensuring monotonic descent. After each accepted swap, we refit Θtm under the updated ordering for consistency. This alternating scheme converges to a locally optimal configuration in practice and typically stabilizes within a few iterations.

- Detailed Description of Matrix Approximation Methods 2D Gaussian. Parameters Θ = {cx , cy , sx , sy , A, b}, where (cx , cy ) is the center, (sx , sy ) are the spreads, A ≥ 0 is the amplitude, and b is an offset. Widths are clipped to a reasonable range to avoid degeneracy. This very compact model captures smooth, bell-shaped structure and aggressively denoises. • 2D Polynomial. Parameters

•

2 Θ = {β0 , β1x , β1y , β2x , β2y , βxy , β3x , β3y , β2xy , βxy }, (48) including an intercept β0 , first-order slopes (β1x , β1y ), second-order curvature terms (β2x , β2y , βxy ), and higher2 order cross-interaction coefficients (β3x , β3y , β2xy , βxy ). 2 We use a low-degree bivariate basis on [0, 1] (e.g., a Legendre-product basis for numerical stability) and apply ℓ2 regularization to mitigate overfitting. This family captures smooth global trends with a moderate parameter count. • Radial Basis Function (RBF). Parameters

approximation bias for substantial compression. When privacy noise is added, the final approximation error satisfies e ≤ ∥A − Ak ∥ + ∥Ak − P eQ e ⊤∥ , ∥A − A∥ | {z } | {z } truncation

(53)

noise & reconstruction

explicitly showing how k (bias) and noise scale (variance) jointly determine fidelity. Thus, low-rank SVD provides scalability and inherent denoising at the utility cost of a controlled approximation bias. Surrogate Parameterization and Parameter Perturbation We present both predefined parametric surrogates and lowrank surrogates under a single parameterization. For each local b t (after the local permutation in (21)), we distance/cost block A m fit a compact parameter vector/tensor Θm and a corresponding reconstruction map R such that b t ≈ R(Θt ). A m m

(54)

We then apply privacy-preserving perturbation directly in parameter space, rather than perturbing matrix entries. This preserves the surrogate’s analytic structure and helps avoid unwhere αk are basis weights, (µkx , µky ) are centers, σ is desirable artifacts (e.g., negative distances or loss of symmetry) a shared bandwidth, and b is an offset. Centers are fixed that can arise from entrywise noise. on a coarse grid with shared bandwidth; only the weights b m 7→ Θm denote the deterministic parameterLet Φm : A are optimized. This yields a localized, multi-peak approxb estimation map, and let ∆Θ = maxA b m ,A b ′ ∥Φm (Am ) − m imation that remains stable under mild regularization. b ′ )∥1 be its global ℓ1 sensitivity under our adjacency Φm (A m notion. Given privacy budget ϵ > 0, we release - Low-Rank SVD ∆Θ We approximate each coefficient matrix A ∈ Rm×n with a e Θm = PD Θm + η m , ηm,j ∼ Lap 0, , (55) ϵ rank-k surrogate that preserves its dominant structure while drastically reducing dimensionality. Compute the SVD A = where P projects the noisy parameters to a valid domain D UΣV⊤ , where U ∈ Rm×m and V ∈ Rn×n are orthogonal, D (e.g., positivity constraints, bounded coefficients, or nonand Σ = diag(σ1 , . . . , σr ) with σ1 ≥ · · · ≥ σr > 0 and r = negativity constraints), ensuring numerical stability and wellrank(A). The truncated SVD retains the top-k components, posed reconstruction. The server reconstructs a perturbed Θ = {α1 , . . . , αK , µ1x , µ1y , . . . , µKx , µKy , σ, b}, (49)

Ak = Uk Σk Vk⊤ , m×k

k ≪ min{m, n},

coefficient block via

(50)

n×k

where Uk ∈ R , Vk ∈ R , and Σk = diag(σ1 , . . . , σk ). For compact storage and upload, we use the factorized parameterization Q := Vk ∈ Rn×k , (51) so that Ak = PQ⊤ . This reduces the representation from mn entries to k(m+n) parameters and concentrates subsequent privacy noise on a small set of variables. Θ = (P, Q),

P := Uk Σk ∈ Rm×k ,

em = R Θ em . A

(56)

a) Instantiation A: predefined parametric functions.: In this case, Θm = {θm,1 , . . . , θm,d } ∈ Rd parameterizes a predefined family f (·; Θ) (2D Gaussian, 2D Polynomial, or RBF), and the reconstruction map is evaluation on the index grid: Rpara (Θm ) ij ≜ f (i, j; Θm ). (57) The projection PD enforces family-specific validity constraints (e.g., σ > 0, bounded polynomial coefficients, or A ≥ 0) so that the released surrogate remains stable and interpretable. b) Instantiation B: low-rank surrogates (SVD/factorization).: In this case, Θm encodes a rank-k factorization, and R composes the low-rank matrix. Two equivalent parameterizations are convenient: (i) SVD factors: Θm = (Uk , Σk , Vk ) with Rsvd (Θm ) = Uk Σk Vk⊤ .

Error guarantees and rank selection. By the Eckart–Young theorem, Ak is the best rank-k approximation to A in both spectral and Frobenius norms: X ∥A − Ak ∥2F = σi2 , ∥A − Ak ∥2 = σk+1 . (52) i>k

A sets k to capture a target energy fraction Pchoice Pkcommon r 2 2 i=1 σi i=1 σi ≥ τ (e.g., τ = 0.95), trading a controlled

TABLE VII. Comparison of different matrix approximation methods. Surrogate 2D Gaussian

Params in Θtm cx , cy , sx , sy , a, b

2D Polynomial

β0 , β1x , β1y , β2x , 10 β2y , βxy , β3x , β3y , 2 β2xy , βxy α1 , . . . , α K , 25 µ1x , µ1y , . . . , µKx , µKy , σ, b Uk , Σk , Vk 30–50

RBF Low-Rank SVD

Size of Θtm 6

Accuracy Low

Notes Compact, but only fits symmetric bell-shaped data well.

Medium

Simple to fit; more effective for smooth surface trends.

High

Flexible and localized, works well for multi-peak or irregular matrices. Captures global structure with fewer components; efficient to store and reconstruct.

Very High

Noise Strategy Subtract Laplace noise from amplitude (a) and increase spread (σ) to lower the surface everywhere. Subtract Laplace noise from coefficients, especially the intercept, to bias the surface downward. Add negative-biased noise to weights (αi − |Lap(b)|); optionally perturb centers for robustness. Subtract noise from singular values (Sk ) or from the reconstructed matrix directly to ensure conservative output.

(ii) Rank-k factors: Θm = (P, Q) with P := Uk Σk ∈ coefficient block in the permuted index space used during local Rm×k and Q := Vk ∈ Rn×k , so that Rfac (Θm ) = fitting. The reconstruction follows the corresponding surrogate t t e e PQ⊤ . family: Am = f ·; Θm , t ∈ T . Because k(m+n) ≪ mn, injecting noise into Θm substantially Instantiation A: predefined parametric reconstruction. reduces posterior leakage relative to entrywise perturbation When a predefined parametric family is used, reconstruction is while preserving approximation fidelity. performed element-wise by evaluating the recovered surface: c) A practical perturbation rule for SVD surrogates.: For e e Rpara (Θ) ≜ f (i, j; Θ). (59) distance and utility-loss blocks, we perturb only the singular ij values while keeping Uk and Vk fixed to maximize utility Concretely, typical instances include per unit of privacy. Let s = diag(Σk ) ∈ Rk≥0 . We draw  2 (y−c̃ )2 x) i.i.d.Laplace noise ηi ∼ Lap(0, b) and apply a downward − 2s̃2y + b̃, Gaussian, Ã exp − (x−c̃  2 2s̃  x y  biased soft-thresholding:   X β̃pq Lp (x)Lq (y) + b̃, Polynomial, e = s̃i = max{si − τ − |ηi |, 0}, i = 1, . . . , k. (58) f (x, y; Θ)  p,q X    Here b = ∆ϵs with ∆s a Lipschitz bound (metric sensitivity) α̃k ϕ (∥[x, y] − µ̃k ∥; σ̃) + b̃, RBF,   of the singular-value map under our adjacency notion; when k (60) a tight analytical ∆s is unavailable, we select b via a binary where L (·) denotes Legendre polynomials (chosen for numerp search in our posterior-leakage auditor to satisfy the target (ϵ, δ)-PmDP constraint. The deterministic shrinkage τ = cϵ ical stability and orthogonality), and ϕ(·; σ̃) is the radial basis b t is then mapped back to the original (with c ∈ [0.3, 1.0] tuned on a small validation set) reduces the kernel. The resulting A m ˜ index order via Eq. (36). probability of overestimation (e.g., d > d) that would otherwise relax mDP constraints, thereby lowering violation rates without Instantiation B: low-rank reconstruction. When a low-rank additional privacy spend. surrogate is used, the coordinator receives perturbed rank-k e k = diag(s̃) and truncate to the conforming factors Θ e = (P, e Q) e (or equivalently (U, e Σ, e V)) e and forms We then form Σ e = Uk (:, 1 : r), V e = the raw reconstruction rank r = #{i : s̃i > 0}, yielding U e = diag(s̃1:r ). The uploaded object is Vk (:, 1 : r), and Σ b =P eQ e ⊤ (equivalently A b =U eΣ eV e ⊤ ), A (61) e e e e e Q)), e Θm = (U, Σ, V) (or its equivalent (P, and the server⊤ e e e e side reconstruction follows from (56): Am = UΣV (or followed by unpermutation using Eq. (36). em = P eQ e ⊤ ). A Structure-preserving tightening. Because noise may violate algebraic or geometric invariants required by the LP constraints, Coefficient Reconstruction The server-side after perturbation, depending on the surro- the coordinator applies a structure-preserving tightening step. gate family used during fitting, the server reconstructs each This yields more conservative (constraint-safe) coefficients and coefficient block through a predefined parametric mapping or empirically reduces the violation ratio. • Metric closure for intra-neighborhood distance blocks. through a low-rank factorization and then applies a structuree dx , we compute the metric closure (Floyd–Warshall): For A m preserving post-processing step (which is privacy-free due to the post-processing immunity) to ensure feasibility and improve e dx (i, j) ← min A e dx (i, k) + A e dx (k, j) . A (62) m m m downstream stability. k∈Nm Server-side reconstruction (parametric and low-rank surThis enforces non-negativity, symmetry, zero diagonal, rogates). The coordinator receives the privatized surrogate and all triangle inequalities by tightening noise-induced e t }t∈T from each user m and reconstructs each parameters {Θ overestimates to the smallest path-consistent upper bounds. m

•

ˆ i , xj ) > d(xi , xj ), then by definition of ∆+ , If d(x

Rectangular triangle tightening for cross blocks. For e m ∈ {A e dy e cy each rectangular block B m , Am }, we tighten entries using triangle-inequality-consistent upper bounds against the metric core: e m (i, j) ← min A e dx (i, k) + B e m (k, j) . B (63) m

ˆ i , xj ) − d(xi , xj ) ≤ ∆+ ≤ γ. d(x Thus, ˆ i , xj ) − γ ≤ d(xi , xj ). d(x

k∈Nm

(66)

Since d¯γ is obtained by clipping dˆ − γ at zero, the inequality still holds. Hence d¯γ ≤ d elementwise. The mDP constraint has the form zi,k ≤ eϵd(xi ,xj ) zj,k .

e dx into the This propagates the conservative geometry of A m rectangular blocks, shrinking overestimates and reducing violations without additional privacy expenditure. e t }t∈T are then used in The tightened reconstructions {A m the global LP exactly where the original coefficient blocks {Atm }t∈T would appear.

Because d¯γ (xi , xj ) ≤ d(xi , xj ), we have ¯

eϵdγ (xi ,xj ) ≤ eϵd(xi ,xj ) .

D. Additional Experimental Results Open Science We provide a complete MATLAB implementation of the proposed DistA-LP framework to enable full reproducibility of all experimental results reported in the paper. The artifact is available at https://anonymous.4open.science/r/DistA-LP-D2F0. It includes a README that documents the code structure, key components, and all relevant configuration parameters.

(67)

Therefore, constraints enforced using d¯γ are at least as strict as those under the true distance d, giving a conservative mDP guarantee. Why local permutation helps.

Let Πm be a user-private random permutation over the neighborhood indices, applied to coefficient blocks as in (15). From the server’s view, the uploaded surrogate parameters are drawn from ah permutation-mixture distribution Pr[Θ̃m ∈ i Violation Ratio S | Am ] = EΠm Pr[Θ̃m ∈ S | Πm (Am )] . In particular, any Table III summarizes the empirical mDP violation ratios for two coefficient tensors that differ only by a relabeling of all perturbation mechanisms across datasets, privacy budgets, neighborhood indices (i.e., A′m = σ(Am ) for some permutation and domain sizes. Predefined-noise baselines (EM and Laplace) σ) induce the same distribution over uploads, since Πm ◦ σ and the hybrid method RMP generally exhibit zero or near- is identically distributed to Πm . Thus, permutation primarily zero violations. In contrast, centralized LP consistently incurs obscures the alignment between indices and concrete locations, the largest violation ratios among all methods, as its coarse- weakening linkage-based inference. grained distance approximation can substantially mis-estimate inter-record distances and thereby relax the enforced mDP constraints. DistA-LP, on the other hand, achieves empirically negligible violations, highlighting that privacy-aware coefficient approximation combined with probabilistic auditing can maintain near-zero mDP violations even in regimes where centralized LP accumulates substantial violations across scales and budgets. Conservative distance repair. Let d denote the true metric and let dˆ denote the reconstructed distance obtained from the perturbed surrogate parameters. For a repair parameter γ ≥ 0, define ˆ i , xj ) − γ, 0}. d¯γ (xi , xj ) = max{d(x (64) Let ˆ i , xj ) − d(xi , xj ))+ . ∆+ = max(d(x i,j

If γ ≥ ∆+ , then d¯γ (xi , xj ) ≤ d(xi , xj ),

∀i, j.

Therefore, any perturbation matrix satisfying the mDP constraints under d¯γ also satisfies the corresponding mDP constraints under the true metric d. ˆ i , xj ) ≤ d(xi , xj ), then Proof. For any pair (i, j), if d(x ˆ i , xj ) ≤ d(xi , xj ). d¯γ (xi , xj ) ≤ d(x

(65)

TABLE VIII. Performance of coefficient approximation across the three datasets (4000 nodes). Distance matrix Ãdx m RelErr ViolRatio

Distance matrix Ãd Nm ×Y RelErr ViolRatio

Utility loss matrix ÃcNm ×Y PL RelErr

Methods

2D Gaussian 2D Polynomial RBF Low-Rank SVD

0.65±0.19 4.02±0.49 7.01±1.66 0.80±0.15

0.56±0.08 380.54±84.93 0.57±0.05 0.07±0.00

0.38±0.08 0.55±0.06 0.41±0.06 0.13±0.02

Rome road map 0.37±0.08 2.48±0.77 7.00±3.59 0.53±0.20

0.55±0.06 213.81±67.31 0.58±0.08 0.05±0.01

0.42±0.08 0.47±0.10 0.39±0.05 0.20±0.02

0.59±0.29 2.30±0.92 6.04±2.97 0.62±0.11

0.60±0.06 48.45±35.55 0.71±0.03 0.03±0.00

2D Gaussian 2D Polynomial RBF Low-Rank SVD

1.29±0.65 2.80±0.79x 3.79±1.04 0.84±0.08

0.56±0.06 338.24±310.64 0.53±0.02 0.07±0.00

London road map 0.37±0.06 0.49±0.15 0.43±0.18 2.89±1.15 0.43±0.02 3.49±0.80 0.14±0.01 0.58±0.04

0.55±0.12 270.3±99.21 0.52±0.02 0.04±0.00

0.44±0.09 0.46±0.09 0.35±0.04 0.21±0.01

0.43±0.14 2.73±0.65 3.40±0.84 0.96±0.39

0.51±0.04 120.87±30.61 0.62±0.01 0.03±0.00

2D Gaussian 2D Polynomial RBF Low-Rank SVD

0.48±0.14 2.22±0.69 2.52±0.77 0.95±0.18

0.66±0.11 268.42±160.48 0.65±0.05 0.08±0.00

New York City road map 0.40±0.09 0.55±0.15 0.67±0.03 0.43±0.21 2.20±1.33 271.57±158.07 0.44±0.03 3.08±1.00 0.64±0.08 0.08±0.00 0.65±0.13 0.07±0.01

0.46±0.06 0.41±0.18 0.45±0.01 0.14±0.01

0.39±0.12 2.28±0.80 2.93±0.89 0.86±0.25

0.54±0.03 63.64±58.32 0.66±0.01 0.03±0.00

TABLE IX. Performance of coefficient approximation across the three datasets (6000 nodes). Distance matrix Ãdx m RelErr ViolRatio

Distance matrix Ãd Nm ×Y RelErr ViolRatio

Utility loss matrix ÃcNm ×Y PL RelErr

Methods

2D Gaussian 2D Polynomial RBF Low-Rank SVD

1.17±1.33 3.43±1.05 3.62±1.45 0.79±0.16

0.60±0.12 350.98±217.18 0.58±0.01 0.07±0.00

0.39±0.08 0.50±0.14 0.41±0.03 0.13±0.01

Rome road map 0.83±0.62 2.58±1.37 5.55±2.90 0.65±0.32

0.50±0.07 207.81±143.53 0.53±0.03 0.04±0.00

0.47±0.07 0.42±0.04 0.36±0.05 0.19±0.02

0.85±0.76 3.02±1.87 3.54±1.19 0.78±0.23

0.56±0.04 6.34±10.45 0.69±0.02 0.04±0.00

2D Gaussian 2D Polynomial RBF Low-Rank SVD

0.55±0.12 2.68±1.25 4.45±2.92 0.87±0.06

0.57±0.04 393.66±152.71 0.53±0.02 0.07±0.00

London road map 0.39±0.05 0.72±0.22 0.51±0.16 2.17±0.71 0.41±0.05 4.45±1.72 0.14±0.01 0.56±0.06

1.17±1.21 463.08±279.57 0.52±0.01 0.04±0.00

0.47±0.06 0.46±0.11 0.39±0.03 0.22±0.02

0.68±0.47 2.97±1.67 3.69±1.87 0.84±0.46

0.45±0.04 85.62±45.48 0.62±0.01 0.03±0.00

2D Gaussian 2D Polynomial RBF Low-Rank SVD

0.54±0.12 2.35±0.59 2.44±0.67 0.96±0.23

0.68±0.06 483.97±284.02 0.66±0.03 0.08±0.01

New York City road map 0.44±0.05 0.72±0.33 1.74±2.24 0.50±0.15 2.56±1.40 352.18±168.12 0.47±0.05 3.30±1.71 0.65±0.04 0.09±0.02 0.65±0.16 0.06±0.01

0.50±0.03 0.43±0.09 0.48±0.08 0.15±0.01

0.65±0.36 2.39±0.93 3.00±1.37 0.91±0.25

0.51±0.05 12.02±15.26 0.67±0.03 0.03±0.00

TABLE X. Performance of coefficient approximation across the three datasets (8000 nodes). Methods

Low-Rank SVD

0.60±0.06

Low-Rank SVD Low-Rank SVD

0.93±0.07 0.86±0.09

Distance matrix Ãdx m RelErr ViolRatio

Distance matrix Ãd Nm ×Y RelErr ViolRatio

Utility loss matrix ÃcNm ×Y PL RelErr

0.07±0.00

0.14±0.01

Rome road map 0.58±0.15

0.04±0.01

0.21±0.01

1.09±0.35

0.04±0.00

0.07±0.00

London road map 0.14±0.01 0.68±0.11

0.04±0.00

0.23±0.01

0.79±0.26

0.04±0.00

0.08±0.00

New York City road map 0.08±0.01 0.67±0.16 0.06±0.00

0.14±0.03

0.73±0.10

0.03±0.00

TABLE XI. Performance of coefficient approximation across the three datasets (10000 nodes). Distance matrix Ãdx m RelErr ViolRatio

Distance matrix Ãd Nm ×Y RelErr ViolRatio

Utility loss matrix ÃcNm ×Y PL RelErr

Methods

Low-Rank SVD

0.59±0.09

0.07±0.00

0.14±0.01

Rome road map 0.62±0.29

0.04±0.00

0.21±0.01

1.09±0.21

0.04±0.00

Low-Rank SVD

0.89±0.10

0.07±0.00

London road map 0.13±0.01 0.56±0.04

0.04±0.00

0.20±0.01

0.80±0.27

0.04±0.00

0.08±0.01

New York City road map 0.07±0.01 0.54±0.12 0.06±0.01

0.15±0.03

0.81±0.27

0.03±0.00

Low-Rank SVD

0.94±0.12

TABLE XII. Memory cost. Number of records Privacy budget (km

−1

2,000 )

ϵ = 4.0

4,000

ϵ = 7.0

ϵ = 10.0

ϵ = 4.0

6,000

ϵ = 7.0

ϵ = 10.0

ϵ = 4.0

ϵ = 7.0

ϵ = 10.0

1.28 MB 1.92 MB 1.92 MB ————>204.7 TB 392 MB 933.56 MB 6.6318 MB

1.28 MB 1.92 MB 1.92 MB ————>204.7 TB 392 MB 933.56 MB 19.645 MB

1.92 MB 2.88 MB 2.88 MB ————>691.1 TB 392 MB 2313.32 MB 1.6471 MB

1.92 MB 2.88 MB 2.88 MB ————>691.1 TB 392 MB 2313.32 MB 8.4639 MB

1.92 MB 2.88 MB 2.88 MB ————>691.1 TB 392 MB 2313.32 MB 25.628 MB

Pre-defined EM Noise Distribution Laplace Hybrid RMP Method COPT Optimization LP Based LP-A Methods PAnDA † D ISTA-LP

0.64 MB 0.96 MB 0.96 MB ————>25.6 TB 392 MB 232.39 MB 1.0375 MB

0.64 MB 0.96 MB 0.96 MB ————>25.6 TB 392 MB 232.39 MB 7.7995 MB

Rome road map 0.64 MB 1.28 MB 0.96 MB 1.92 MB 0.96 MB 1.92 MB ————————>25.6 TB >204.7 TB 392 MB 392 MB 232.39 MB 933.56 MB 14.58 MB 1.4009 MB

Pre-defined EM Noise Distribution Laplace Hybrid RMP Method COPT Optimization LP Based LP-A Methods PAnDA † D ISTA-LP

0.64 MB 0.96 MB 0.96 MB ————>25.6 TB 522.72 MB 401.42 MB 1.7227 MB

0.64 MB 0.96 MB 0.96 MB ————>25.6 TB 522.72 MB 401.42 MB 7.8758 MB

0.64 MB 0.96 MB 0.96 MB ————>25.6 TB 522.72 MB 401.42 MB 22.355 MB

London road map 1.28 MB 1.28 MB 1.92 MB 1.92 MB 1.92 MB 1.92 MB ————————>204.7 TB >204.7 TB 522.72 MB 522.72 MB 923.15 MB 923.15 MB 2.3711 MB 8.7055 MB

1.28 MB 1.92 MB 1.92 MB ————>204.7 TB 522.72 MB 923.15 MB 22.77 MB

1.92 MB 2.88 MB 2.88 MB ————>691.1 TB 522.72 MB 2309.97 MB 2.7703 MB

1.92 MB 2.88 MB 2.88 MB ————>691.1 TB 522.72 MB 2309.97 MB 11.913 MB

1.92 MB 2.88 MB 2.88 MB ————>691.1 TB 522.72 MB 2309.97 MB 23.616 MB

Pre-defined EM Noise Distribution Laplace Hybrid RMP Method COPT Optimization LP Based LP-A Methods PAnDA D ISTA-LP†

0.64 MB 0.96 MB 0.96 MB ————>25.6 TB 184.95 MB 229.57 MB 1.5049 MB

0.64 MB 0.96 MB 0.96 MB ————>25.6 TB 184.95 MB 229.57 MB 4.0234 MB

New York City road map 0.64 MB 1.28 MB 1.28 MB 0.96 MB 1.92 MB 1.92 MB 0.96 MB 1.92 MB 1.92 MB ————————————>25.6 TB >204.7 TB >204.7 TB 184.95 MB 184.95 MB 184.95 MB 229.57 MB 922.61 MB 922.61 MB 16.244 MB 1.6511 MB 5.335 MB

1.28 MB 1.92 MB 1.92 MB ————>204.7 TB 184.95 MB 922.61 MB 20.137 MB

1.92 MB 2.88 MB 2.88 MB ————>691.1 TB 184.95 MB 2312.96 MB 6.4273 MB

1.92 MB 2.88 MB 2.88 MB ————>691.1 TB 184.95 MB 2312.96 MB 11.796 MB

1.92 MB 2.88 MB 2.88 MB ————>691.1 TB 184.95 MB 2312.96 MB 23.774 MB

TABLE XIII. Communication cost. Surrogate Raw Matrices 2D Gaussian 2D Polynomial RBF Low-Rank SVD

Upload to Server (Bytes)

Download from Server (Bytes)

Upload to Downstream (Bytes)

Total End-to-End (Bytes)

Total End-to-End (kb)

9600 304 400 760 1640

3200 3200 3200 3200 3200

8 8 8 8 8

12808 3512 3608 3968 4848

12.5078 3.4297 3.5234 3.875 4.7344

Runtime Fraction: Swap vs Parameter Fitting

120

Iterations to Stop (Converged or Max Iter)

Average Loss vs Iteration (All Locations)

6600

Best during swaps (swap-only) Refit-only (no swaps) After swap + refit (actual)

6400

100

40 Fit Swap

Average Loss

Count

Percent of total runtime

6200 50

6000

4950

5800 4900

5600 4850

5400

40 20

4800 5

5200

5000

4800

0 0

100

150

200

Location

Outer iteration

Outer iterations executed

(a) Gaussian surrogate Runtime Fraction: Swap vs Parameter Fitting

Iterations to Stop (Converged or Max Iter)

140

105Average Loss vs Iteration (All Locations)

2.5

Best during swaps (swap-only) Refit-only (no swaps) After swap + refit (actual)

120

Average Loss

100

Count

Percent of total runtime

Fit Swap

80 60

5500

1.5

5000 4500 4000

3500 3000

0.5 10

0 0

100

150

200

Location

Outer iteration

Outer iterations executed

(b) Polynomial surrogate Runtime Fraction: Swap vs Parameter Fitting

Iterations to Stop (Converged or Max Iter)

Average Loss vs Iteration (All Locations)

6500

Best during swaps (swap-only) Refit-only (no swaps) After swap + refit (actual)

6000

5500

Average Loss

40 50

Fit Swap

Count

Percent of total runtime

40 30

5000

2400 2350

4500

2300

4000 2250

3500 20

2200 6

3000

10 10

0 0

100

150

2500

200

2000 4

Location

Outer iterations executed

Outer iteration

Fig. 6: Surrogates approximation performance with a maximum of 10 optimization iterations, showing the runtime split between parameter fitting and swap search; convergence/stopping iterations; and the decrease in average loss under swap-only, refit-only, and combined swap–refit updates. Runtime Fraction: Swap vs Parameter Fitting

Iterations to Stop (Converged or Max Iter)

Average Loss vs Iteration (All Locations)

6600

Best during swaps (swap-only) Refit-only (no swaps) After swap + refit (actual)

6400 25

6200 6000

Average Loss

20 Fit Swap

Count

Percent of total runtime

4820

5800

4800 4780

5600

4760

5400

4740 4720

5200 4700 14

5000

4800 0

0 0

100

Location

150

200

Outer iterations executed

(a) Gaussian surrogate

4600 0

Outer iteration

Runtime Fraction: Swap vs Parameter Fitting

Iterations to Stop (Converged or Max Iter)

105

2.5

Average Loss vs Iteration (All Locations) Best during swaps (swap-only) Refit-only (no swaps) After swap + refit (actual)

30 2

25 5500

Fit Swap

Average Loss

Count

Percent of total runtime

1.5

5000 4500 4000

3500

20 10

3000 2500

0.5

0 0

100

150

200

Location

Outer iterations executed

Outer iteration

(b) Polynomial surrogate Runtime Fraction: Swap vs Parameter Fitting

Average Loss vs Iteration (All Locations)

6500

Best during swaps (swap-only) Refit-only (no swaps) After swap + refit (actual)

6000

5500

40 Fit Swap

Average Loss

5000

Count

Percent of total runtime

Iterations to Stop (Converged or Max Iter)

20 15

2400 2350

4500

2300 2250

4000

2200

3500

2150 6

3000

5 2500

0 0

100

150

200

2000

Outer iterations executed

Location

Outer iteration

Fig. 7: Surrogate approximation performance with a maximum of 20 optimization iterations, showing the runtime fraction between the fitting and swap search; convergence/stopping-iteration distribution; and the average loss decay across iterations for swap-only, refit-only, and combined swap–refit updates.

vs mDP Distance-Overestimation Violation Ratio for DLR

0.4

vs mDP Distance-Overestimation Violation Ratio for DLR2obf

0.3

Reconstructed D LR2obf

Reconstructed D LR

Original Violation Ratio

0.35

0.25

0.3

Violation ratio

0.2

0.25

0.2

0.15

0.15 0.1

0.1 0.05

0.05

0 0

0.5

1.5

2.5

3.5

4.5

repair parameter

0.5

1.5

2.5

3.5

4.5

repair parameter

Fig. 8: γ vs mDP distance-overestimation violation ratio for both intra-neighborhood distance block and crossdomain distance block.

A PPENDIX A A RTIFACT A PPENDIX A. Description & Requirements 1) How to access: The artifact is available at https://zenodo.org/records/21782566. It contains the MATLAB implementation, experiment inputs, result directories, and documentation needed to evaluate DISTA-LP and its baselines. The primary entry point is Bendersdecomposition/main.m; main_surrogate.m evaluates coefficient approximation separately, and baseline entry points are under Bendersdecomposition/baseline/. 2) Hardware dependencies: The artifact was evaluated on a workstation with 32 GB RAM. We recommend a multicore CPU, at least 32 GB RAM, and 5–10 GB free storage. A GPU is not required. Begin with a 2,000-node experiment; centralized LP, COPT and larger PAnDA settings can require substantially more memory and time. 3) Software dependencies: The main framework requires MATLAB R2024b or later and Optimization Toolbox. Statistics and Machine Learning Toolbox is required by statistical and sampling functions. Parallel Computing Toolbox is recommended for parallel subproblems. Baselines were checked with MATLAB R2025b; baseline functions using heaviside require Symbolic Math Toolbox. The code supports Windows, macOS, and Linux when MATLAB and the required toolboxes are available. 4) Benchmarks: The supplied experiments use road-network data for Rome, London, and New York City (NYC), with 2,000, 4,000, or 6,000 selected records, privacy budgets ϵ ∈ {4, 7, 10} km−1 , ten users, and dataset-dependent repetitions: Setting

lines 14-17 for the Rome dataset, lines 20-23 for the NYC dataset, and lines 26-29 for the London dataset. b) Target regions: In main.m, activate the Regions block matching the city and record count. The region count R = numel(Regions) also determines the number of repetitions and columns in the uploaded-budget matrix. From lines 124-163, based on the dataset and node count. c) Node count (nr loc selected) and uploaded budget: Set, for example, env parameters.nr loc selected = 2000; budget assigned = budget assigned rome svd. The selected budget assigned must match the city, node count, and surrogate. The main runner uses epsilon value = [10, 7, 4]; so matrix rows 1–3 correspond to ϵ = 10, 7, 4, respectively. With approx matrix=1, the budget remaining for optimization is env parameters.EPSILON = EPSILON - budget assigned(ep idx,r); Setting approx matrix=0 assigns the full ϵ to optimization and disables approximated-matrix privacy accounting. The uploaded budget is between lines 221-263, budget assigned at line 363, and node count at line 383. d) Candidate outputs: In classes/User/User.m, set the fraction used by obf identifier at line 201:

Repetitions

Rome, 2,000 records Rome, 4,000/6,000 records London, all sizes NYC, all sizes

5 4 6 4

Records

Fraction

Count

2,000 4,000 6,000

0.01 0.005 0.0033

20 20 ≈ 20

The implementation retains the nearest candidates by Haversine distance when more locations fall inside the obfuscation range. e) Consistency check.: Confirm that the dataset loader, Regions, nr loc selected, uploaded-budget definition and assignment, and candidate-output fraction all describe the same experiment. Use names and comments rather than relying only on line numbers, which may change. C. Experiment Workflow

The evaluated baselines are Exponential Mechanism (EM), Bayesian remapping (RMP), coarse-grid LP approximation (LPA), Planar Laplace, PAnDA, COPT and the original centralized metric-DP LP. The full centralized LP and COPT are infeasible at the evaluated domain sizes and are reported as unavailable. B. Artifact Installation & Configuration Download the zip folder and enter the framework directory (Bendersdecomposition folder): Start MATLAB in this directory. The main runner adds the required class and function directories to the MATLAB path. Before a main-framework run, select a consistent dataset, target-region block, node count, uploaded-budget matrix, budget assigned, and obfuscated-location fraction. a) Dataset: In the Bendersdecomposition folder, func/read_files/readCityMapInfo.m, leave exactly one of the dataset (Rome, London, or NYC) loading blocks active depending on the dataset of consideration. From

The main workflow loads a road network and sampled user inputs, constructs private surrogate representations of LP coefficients, accounts for privacy budget consumed during parameter upload, solves the distributed problem through Benders decomposition, and aggregates the output metrics. Run it from Bendersdecomposition/: r u n ( ’ main . m’ ) It reports utility loss, mDP violation ratio, computation time, privacy-loss statistics for uploaded parameters, and surrogate approximation errors. The baseline workflow is independent. From Bendersdecomposition/baseline/, edit the configuration near the top of each runner: city = ’london’; node count = 2000; epsilons = [4, 7, 10]; user ids = 1:10; Then run: r u n ( ’EM EMBR. m’ ) % EM, RMP, LP−A r u n ( ’ l a p l a c e b a s e l i n e . m’ ) r u n ( ’ PAnDA baseline . m’ ) % PAnDA and LB r u n ( ’ C O P T b a s e l i n e . m’ )

r u n ( ’ L P b a s e l i n e . m’ ) The baseline runners print tables to the MATLAB Command Window and do not create new MAT or CSV files. Their displayed epsilon order is [4, 7, 10], unlike the main runner’s internal budget-row order [10, 7, 4]. D. Major Claims (C1) Privacy: DISTA-LP enforces the target metric differential privacy constraints with near-zero empirical violation ratios across the evaluated cities and privacy budgets. Validate with E1 and compare with E2. • (C2) Utility: DISTA-LP improves utility over predefined mechanisms and the PAnDA hybrid while remaining above the lower-bound solution. Validate with E1 and E2. • (C3) Scalability: Benders decomposition makes finegrained road-network experiments tractable where full centralized LP/COPT runs are unavailable or exceed the 1,800-second limit. Validate with E1 and E2 across node counts. • (C4) Private coefficient upload: Compact surrogates and calibrated parameter perturbation allow coefficient information to be uploaded with measurable privacy loss and approximation error. Validate with E1.

•

E. Evaluation 1) Experiment (E1): DISTA-LP main framework: Estimated effort: 10–20 human-minutes for setup; compute time is machine- and configuration-dependent. Use 2,000 records for the initial functional test. Preparation. Select one city in readCityMapInfo.m; activate its matching Regions, uploaded-budget matrix, and budget-assigned in main.m; set nr loc selected; set approx matrix=1; and set the matching fraction in User.m. For a Rome 2,000-record test, use five regions, budget assigned rome svd, and fraction 0.01. Execution. From Bendersdecomposition/, run main.m. Allow all active regions and epsilon values to finish. Repeat for each city and record count required by the paper result being checked. Results. Record utility loss, violation ratio, computation time, upload privacy loss, and surrogate relative/violation errors. Interpret uploaded-budget rows in [10, 7, 4] order. Confirm successful optimization and inspect any recorded error message before using a result. E1 evaluates C1–C4. Utility loss is divided by 10,000 before printing. Thus 3 4 7 0 8 . 0 2 +/ − 2 4 1 4 . 5 9 becomes 3 . 4 7 0 8 0 2 +/ − 0 . 2 4 1 4 5 9 , o r 3 . 4 7 +/ − 0 . 2 4 a f t e r r o u n d i n g . 2) Experiment (E2): Baselines: Estimated effort: 10 human-minutes for configuration; compute time varies widely. LP-A, PAnDA, COPT, and centralized LP can be long-running or resource-limited at larger sizes. Preparation. From the baseline directory, set city, node count, epsilons=[4,7,10], and user ids=1:10. Both

COPT_baseline.m and LP_baseline.m use Optimization Toolbox’s linprog, set max time seconds=1800. A oneuser run may be used only as a functional check; restore ten users for paper reproduction. Execution. Run EM_EMBR.m, laplace_baseline.m, PAnDA_baseline.m, COPT_baseline.m and LP_baseline.m. Copy the printed utility-loss, violationratio, and computation-time tables. Repeat for the remaining city/size combinations. Results. For each repetition, the scripts sum utility loss across ten users and average runtime and violation ratio across users. Then report the mean and sample standard deviation across repetitions. EM, RMP, and Planar Laplace display zero violation by construction/post-processing; LP-A reports an empirical distance-approximation violation ratio. PAnDA prints separate PAnDA and companion lower-bound rows. E2 evaluates C1–C3. 3) Experiment (E3): Surrogate approximation.: Estimated effort: 10–20 human-minutes per city; compute time depends on configuration. This experiment reproduces the RelErr and ViolRatio columns of Table II for Gaussian, polynomial, RBF, and low-rank SVD surrogates. Preparation. In func/read_files/readCityMapInfo.m, activate exactly one dataset-loading block: Rome (lines 14–17), NYC (20–23), or London (26–29). In classes/User/User.m, uncomment fitting calls at lines 113–115, noisy-parameter calls at 118–120, and reconstruction calls at 123–125; keep the corresponding SVD calls enabled. Set the obf identifier fraction at approximately line 201 to 0.01. In main_surrogate.m, activate the corresponding 2,000-node region block: Rome (124–128), NYC (143–146), or London (156–161). Ensure the effective nr loc selected is 2000, with NR LOC=1 if the later NR LOC*2000 assignment remains active. Execution. From Bendersdecomposition/, run main_surrogate.m. Allow all regions and users to finish. Preserve Summary_Gaussian.csv, Summary_Polynomial.csv, Summary_RBF.csv, and Summary_SVD.csv before repeating for the other cities. Results. Compare generated summaries with those in rome_relErr_violRatio/, nyc_relErr_ violRatio/, and london_relErr_violRatio/, relative to Bendersdecomposition/. Match Rel_DLR/Viol_DLR, Rel_DLR2/Viol_DLR2, and Rel_CRL to the coefficient blocks in Table II. Surrogate ViolRatio differs from the final mechanism’s metric-DP violation ratio. Compare original and reconstructed coefficients before User.apply svd() overwrites the working matrices; otherwise artificial zero errors can result. The dedicated surrogate runner avoids this replacement. A modified optimization runner must retain original matrices or use the preserved *_test copies. E3 evaluates C4. 4) Experiment (E4): End-to-end result check: Estimated effort: 15 human-minutes after experiments complete. Preparation. Create one table per city and node count with columns ϵ = 4, 7, 10 and rows for DISTA-LP and every

available baseline. Execution. Transfer E1 and E2 values without changing their epsilon ordering. Compare utility loss, violation ratio, and runtime with the corresponding paper table or figure. Results. Reproduced results should preserve the paper’s qualitative ranking. As a supplied reference point, the main README reports Rome, 4,000 records, ϵ = 4: DISTA-LP utility loss 2.88 ± 0.21 in 10,000-meter units, compared with approximately 7.16 ± 2.53 for EM/Laplace, 6.76 ± 2.55 for LP, and a lower bound near 1.7. Treat this value only as a cross-check; the paper and generated outputs are authoritative. F. Customization The artifact can be customized by selecting another supported city, record count, epsilon, user subset, or surrogate family. Root-level scripts test Gaussian, polynomial, RBF, and lowrank SVD surrogates. When adding a new city or size, supply compatible road-network CSVs, sampled MAT inputs, target-region definitions, graph objects where required by baselines, and a budget matrix with one column per repetition. Preserve the epsilon-to-row mapping, document all parameter changes and use separate output directories to avoid overwriting reference runs. G. Notes Baseline output is currently mean ± sample standard deviation, not a 95% confidence interval. If the paper uses √ 95% normal-approximation intervals, compute 1.96s/ R and label the statistic consistently. • Divide both the utility-loss mean and its deviation by 10,000. • Primary baseline runners do not read or overwrite historical MAT files. Auxiliary files containing memory can write profiling checkpoints. Main and surrogate runnerscan overwrite summaries or sampled inputs; preserve eachrun before reconfiguration • main_BD.m is legacy/development code and is not a primary baseline entry point. •

DISTA-LP: Privacy-Preserving Distributed Linear Programming for Metric Differential Privacy

Analyze & explore related research in Conceptio

Related documents