arXiv:2605.10015v1 [stat.ML] 11 May 2026
Differentially Private Sampling from Distributions via Wasserstein Projection
Shokichi Takakura LY Corporation [email protected]
Seng Pei Liew LY Corporation [email protected]
Satoshi Hasegawa LY Corporation [email protected]
Abstract In this paper, we study the problem of sampling from a distribution under the constraint of differential privacy (DP). Prior works measure the utility of DP sampling with density ratio-based measures such as KL divergence. However, such formulations suffer from two key limitations: 1) they fail to capture the geometric structure of the support, and 2) they are not applicable when the supports of the distributions differ. To deal with these issues, we develop a novel framework for DP sampling with Wasserstein distance as the utility measure. In this formulation, we propose Wasserstein Projection Mechanism (WPM), a minimax optimal mechanism based on Wasserstein projection. Furthermore, we develop efficient algorithms for computing the proposed mechanisms approximately and provide convergence guarantees.
1
Introduction
Differential privacy (DP) [Dwork et al., 2006, Dwork and Roth, 2014] provides a rigorous framework for quantifying and ensuring privacy in data analysis and machine learning. Especially, local differential privacy (LDP) [Kasiviswanathan et al., 2011] is a strong notion of privacy that does not require a trusted server and thus has been widely adopted in practice [Erlingsson et al., 2014, Differential Privacy Team, 2017]. However, most of the prior work on LDP assumes that each user holds only a single data point, which is unrealistic since modern devices often collect and store a large amount of data [Husain et al., 2020]. Several works [Levy et al., 2021, Ghazi et al., 2021] consider the user-level privacy by assuming that each user holds a dataset of the same size generated from an underlying distribution. Although this scenario is more realistic, the requirement of equal dataset size across users is still restrictive. To mitigate this issue, a line of work [Husain et al., 2020, Park et al., 2024, Zamanlooy et al., 2025, Ghoukasian et al., 2025] has treated each user’s dataset as a probability distribution and proposed differentially private mechanisms, which take a distribution as input and release a single sample. We refer to this problem as private sampling. Some of these works have derived optimal mechanisms in terms of worst-case f -divergence between the input and output distributions for various settings [Park et al., 2024, Zamanlooy et al., 2025, Ghoukasian et al., 2025]. While these works have provided valuable insights into the design of private sampling mechanisms, most of them rely on density ratio-based measures such as KL divergence to quantify the utility of the mechanism. However, such formulations suffer from two key limitations: Preprint.
• Support Mismatch: Metrics based on density ratios are not meaningful when the supports of the distributions differ. • Geometric Ignorance: They do not capture the geometric structure of the underlying space (support), which can be crucial for downstream applications. For example, location data is often represented as points in a metric space [Chatzikokolakis et al., 2015], and in recommendation systems and natural language processing, items and words are typically embedded in a latent space where distances capture semantic similarity [Koren et al., 2009, Mikolov et al., 2013], making the geometry of the space crucial for evaluating utility. Moreover, support mismatch arises naturally in these settings: geographic locations may be aggregated to representative points [Chatzikokolakis et al., 2015], a recommender may sample from a broader catalog beyond observed items [Ricci et al., 2010], and language models often restrict or modify the candidate set (e.g., via truncation or special tokens) [Vinod et al., 2025]. To address these limitations, we propose a novel framework for private sampling that optimizes the worst-case Wasserstein distance [Villani et al., 2009] between the privatized distribution of the mechanism and the target distribution. Since Wasserstein distance captures the geometric structure of the support and is properly defined even when the supports of the distributions differ, our framework provides a more flexible way to design private sampling mechanisms. Furthermore, we develop a Wasserstein Projection Mechanism (WPM) based on Wasserstein projection, which is shown to be minimax optimal in terms of worst-case Wasserstein distance. Since computing the exact projection can be computationally expensive, we propose a practical algorithm for computing the approximate projection using entropic regularization and alternating projections. We provide linear convergence guarantees and approximation error bounds for our algorithm. Finally, we demonstrate the effectiveness of our proposed method through numerical experiments on synthetic and real-world datasets. Our contributions are as follows: • To mitigate the limitations of density ratio-based formulations, we propose a novel framework for designing private sampling mechanisms by utilizing Wasserstein distance as the utility measure. • We develop Wasserstein Projection Mechanism (WPM), a minimax optimal mechanism in terms of worst-case Wasserstein distance. • We propose practical algorithms for computing the optimal mechanism using entropic regularized Wasserstein distance and alternating projections. By exploiting the special structure of the problem, we derive a linear convergence guarantee for our algorithm. 1.1
Related Work
Statistics Estimation in LDP A rich literature studies estimating statistics from data privatized under local differential privacy [Asi et al., 2022], including frequency/histogram and mean estimation. Foundationally, Duchi et al. [2013] derived minimax rates for several estimation tasks under LDP, and Kairouz et al. [2014] characterized extremal LDP mechanisms. In contrast, our focus is LDP sampling mechanisms that control distributional distortion (here, in Wasserstein distance) rather than estimating a fixed statistic. This can be viewed as a generalization of the classical estimation problem, where the input distribution is a Dirac measure at a single point. See Appendix A for a more detailed discussion. Wasserstein Distance in DP Several works have utilized the notion of Wasserstein distance in the context of differential privacy [Feldman et al., 2024, Rodríguez-Vítores et al., 2025, Gu et al., 2025]. For instance, Lê Tien et al. [2019], Rakotomamonjy and Liva [2021] proposed randomized mechanisms to compute the (sliced) Wasserstein distance privately. Feldman et al. [2024] utilized Wasserstein distance as a utility measure for private density estimation. In contrast to our work, which considers privacy leakage of a sample from the distribution, in their scenario, the estimated distribution itself is the output of the mechanism. Therefore, while these works have utilized similar terminology, their focus is completely different from ours. 2
1.2
Notation
Let R+ (R++ ) denote the set of non-negative (positive) real numbers, and [k] = {1, 2, . . . , k} for k ∈ N. For x ∈ X , we denote by δx the Dirac measure at x, i.e., δx (A) = 1 if x ∈ A and δx (A) = 0 otherwise. For a ∈ R, we denote by (a)+ the positive part of a, i.e., (a)+ = max{a, 0}.
2
Preliminaries
2.1
Differentially Private Sampling
Differential privacy is a mathematically rigorous framework for quantifying the privacy guarantees of algorithms that operate on sensitive data. Let X be a measurable space equipped with a σ-algebra F, and let P(X ) be the set of all probability measures on X . In the local model of differential privacy (LDP), a client has an input distribution µ ∈ P(X ) which depends on sensitive data. Therefore, sampling directly from µ can potentially leak sensitive information about the data. To ensure privacy, a client samples from a transformed distribution M[µ], where M is a functional mapping P(X ) to P(X ). We require that M satisfy the constraints of differential privacy, which ensures that a sample from the distribution M[µ] does not reveal too much information about the input distribution µ while preserving as much utility as possible. Formally, we have the following definitions of LDP. Definition 2.1 (Local Differential Privacy (LDP)). A mechanism M is said to satisfy ε-differential privacy for some ε > 0, if for all pairs of input distributions µ, µ′ ∈ P(X ) and for all measurable subsets S ⊆ X , we have: M[µ](S) ≤ eε M[µ′ ](S). Conventionally, ε-DP is defined for randomized mechanisms and here, we identify a functional M : P(X ) → P(X ) with a randomized mechanism that samples from M[µ]. 2.2
KL Projection Mechanism
Here, we briefly review the f -divergence based mechanism proposed by Husain et al. [2020], Park et al. [2024], which we refer to as KL projection mechanism (KPM), for private sampling under LDP. For a convex function f : R+ → R with f (1) = 0, the f -divergence between two probability measures µ and ν is defined as Z dµ Df (µ∥ν) = f dν, dν X where dµ dν is the Radon-Nikodym derivative of µ with respect to ν. Common examples of f divergences include the Kullback-Leibler (KL) divergence and χ2 -divergence. In the discrete setting X = [k] for some k ∈ N, the KL projection mechanism [Park et al., 2024] is defined as follows: 1 1 KL µ(x), ε , M [µ](x) = max rµ e +k−1 where rµ > 0 is a normalization constant so that MKL [µ] is a probability distribution. This is minimax optimal in terms of worst-case f -divergence supµ Df (µ∥M[µ]) for any f -divergence [Park et al., 2024]. Furthermore, this can be regarded as a projection of µ with respect to KL divergence ε/2 onto the LDP polytope with base measure m(x) = eεe+k−1 defined as follows: Definition 2.2 (LDP Polytope). For a privacy parameter ε, the LDP polytope Qm,ε is defined as Qm,ε = {ν ∈ P(X ) : e−ε/2 m(A) ≤ ν(A) ≤ eε/2 m(A) for all measurable A ⊆ X },
(1)
where m is a base measure on X . Note that m is not necessarily a probability measure. In the following, we only consider base measures m such that Qm,ε ̸= ∅, i.e., e−ε/2 m(X ) ≤ 1 ≤ eε/2 m(X ). It is not straightforward to extend this mechanism to continuous settings because f -divergences can be infinite when µ is not absolutely continuous with respect to ν. Prior works deal with this issue by restricting the input space of the mechanism. See Park et al. [2024] for details. 3
While f -divergences are widely used in the literature to measure the distance between probability distributions, they have several limitations. First, f -divergences are not meaningful when the supports of the distributions differ. For instance, in the case where we can only sample from a restricted output space V ⊂ X , µ may not be absolutely continuous with respect to M[µ] and thus, the Radondµ Nikodym derivative dM[µ] does not exist 1 . Second, they do not capture the geometric structure of the underlying space X . Let us consider the case where the input distribution µ is a Dirac measure at 1 a point x ∈ X , i.e., µ = δx . In this case, KPM assigns a positive mass eε +k−1 to all points except x including points far from x in a certain metric, which can lead to a large distortion.
3
Proposed Method: Wasserstein Projection Mechanism
To address these limitations of KPM, we propose a novel framework for designing private sampling mechanisms by utilizing Wasserstein distance as the utility measure. Then, we develop a Wasserstein Projection Mechanism (WPM) based on Wasserstein projection, which is shown to be minimax optimal in terms of worst-case Wasserstein distance. 3.1
Wasserstein Distance
In this paper, we assume that the underlying space X is equipped with a metric d : X × X → R+ . For instance, if X is a subset of Rd , we can use the Euclidean distance as the metric. In such a case, the Wasserstein distance is a natural choice for measuring the distance between probability distributions on X , as it captures the geometric structure of the space. For any µ ∈ P(X ) and ν ∈ P(X ), the Wasserstein distance of order p is defined as Wp (ν, µ) =
Z
d(x, y)p dπ(x, y)
inf π∈Π(ν,µ)
1/p ,
X ×X
where Π(ν, µ) is the set of all couplings of ν and µ, i.e., the set of all probability measures on X × X with marginals µ and ν. For simplicity, we assume that X is a compact Polish space, which ensures that the Wasserstein distance is properly defined and finite for all probability measures supported on X . This is in contrast to density ratio-based divergences (e.g., KL), which can be infinite when the supports differ. 3.2
Generalized Problem Formulation
In this paper, we generalize the problem formulation of prior works by allowing the supports of the input and output distributions to differ. Specifically, let V ⊆ X be a measurable subset of X representing the output space of the mechanism M. Then, we aim to design a mechanism M : P(X ) → P(V) that satisfies ε-LDP and minimizes the Wasserstein distance between M[µ] and µ. In particular, we consider the following worst-case cost: U [M] =
sup Wp (M[µ], µ). µ∈P(X )
This general formulation allows us to consider a sampling mechanism from a restricted output space V, which can be useful in practice when we want to limit the range of the output distribution for computational or privacy reasons. Note that this general formulation cannot be handled with density ratio-based formulations, since f -divergences do not provide a meaningful distance when the supports of the distributions differ. 3.3
Wasserstein Projection Mechanism
To solve the above problem, we propose a Wasserstein Projection Mechanism (WPM) defined as follows: 1 Total variation distance can be defined even when the supports differ but discrepancies outside the common support mainly
capture support mismatch and are not substantively meaningful.
4
Definition 3.1 (Wasserstein Projection Mechanism). Wasserstein Projection Mechanism with base measure m and privacy parameter ε is defined as Mm,ε [µ] ∈ arg min Wp (ν, µ), ν∈Qm,ε
where Qm,ε is the LDP polytope defined as in Eq. (1). Since Qm,ε is convex and closed, a minimizer exists under our assumptions. This mechanism satisfies ε-LDP as shown in the following proposition: Proposition 3.2. The Wasserstein Projection Mechanism Mm,ε satisfies ε-LDP. See Appendix C for the proof. Furthermore, we can show that given any ε-LDP mechanism M, there exists a Wasserstein Projection Mechanism Mm,ε that has a uniformly better or equal utility across all input distributions µ. Proposition 3.3 (Uniform Optimality of WPM). For any ε-LDP mechanism M, there exists a base measure m such that for all µ ∈ P(X ), Wp (M[µ], µ) ≥ Wp (Mm,ε [µ], µ). See Appendix D for the proof. The key ingredient of this proof is the fact that M satisfies ε-LDP is equivalent to the range of M being contained in LDP polytope Qm,ε for some base measure m. While prior work [Husain et al., 2020, Park et al., 2024] utilizes the LDP polytope as a sufficient condition for satisfying LDP, we show that this is also a necessary condition. This theorem implies a strong optimality guarantee for the WPM Mm,ε : given any ε-LDP mechanism M, we can construct a WPM Mm,ε that has a uniformly better or equal utility across all input distributions µ. In other words, as long as the Wasserstein distance is used as the utility, there is no point in using other mechanisms. 3.4
Minimax Optimality and Optimal Base Measure
In the following, we focus on the finite discrete case, i.e., X = [k] for some k ∈ N. Without loss of generality, we assume that V = [kv ] for some kv ≤ k. Due to the uniform optimality guarantee of the WPM, WPM with the optimal base measure m∗ is minimax optimal in the sense that it minimizes the worst-case cost across all ε-LDP mechanisms. Theorem 3.4 (Minimax Optimality). There exists a base measure m∗ such that U [M] ≥ U [Mm∗ ,ε ] for any ε-LDP mechanism M. See Appendix E for the proof. The optimal base measure m∗ can be obtained by optimizing the worst-case cost U [Mm,ε ]. In the case of f -divergences, the optimal base measure is uniform over the support of the target distribution since they are agnostic to the geometry of the space X . On the other hand, the optimal base measure for WPM depends on the geometry of the space X , and thus, finding the optimal base measure is a non-trivial problem. Since evaluating the worst-case cost U [M] involves 1) calculating the Wasserstein distance for each input distribution µ and 2) taking the supremum over all input distributions µ, it is computationally intractable to directly solve the above minimax optimization problem to find the optimal base measure. To mitigate this issue, we show that U [Mm,ε ] admits a tractable reformulation and interestingly, it is convex in m, which allows us to efficiently find the optimal base measure via convex optimization techniques. Theorem 3.5 (Simpler Reformulation). Let Cij = d(i, j)p for i ∈ [k], j ∈ [kv ], and α := Pkv e−ε/2 , β := eε/2 . For each m ∈ Rk+v (1/β ≤ j=1 mj ≤ 1/α) and i ∈ [k], define ϕi (m) = τi (m) +
kv X
(α(Cij − τi (m))+ − β(τi (m) − Cij )+ )mj ,
j=1
5
Pkv where τi (m) := min{t : Then, f (m) := j=1 (βmj 1{Cij ≤t} + αmj 1{Cij >t} ) ≥ 1}. maxi∈[k] ϕi (m) is a convex function of m, and for any ε-LDP mechanism M, we have U [M] ≥ f (m∗ )1/p = U [Mm∗ ,ε ], where m∗ ∈ arg minm∈Rkv :1/β≤P mj ≤1/α f (m). +
j
See Appendix F for the proof. This reduction is based on the following key observations: 1) the worst-case utility is achieved by a Dirac measure at some point i ∈ [k], 2) the optimal projection cost v for a Dirac measure has a simple closed-form expression. Once the costs {Cij }kj=1 are sorted for each i, the function ϕi (m) can be evaluated in O(k) time and thus, f (m) can be evaluated in O(k 2 ) time. In addition, f (m) is a convex function of m, and thus, the optimal base measure can be obtained via standard mirror descent algorithms. From a standard argument in convex optimization, we can obtain an εopt -accurate solution in O(log k/ε2opt ) iterations. Algorithmic details and convergence analysis are given in Appendix G. Extension to continuous spaces While we focus on the finite discrete case for computational tractability, our framework is not limited to discrete spaces and can be applied to continuous settings as well. As a concrete example, we derive the optimal base measure when X is the d-dimensional unit sphere S d equipped with the Euclidean distance by exploiting the rotational symmetry of the sphere. See Proposition H.1 in Appendix H for details.
4
Practical Implementation
Given a base measure m and privacy parameter ε, computing the Wasserstein projection Mm,ε [µ] for a given input distribution µ can be formulated as the following Linear Program (LP): min
kv k X X
k×kv v ν∈Rk + , π∈R+ i=1 j=1
Cij πij s.t.
kv X
πij = µi for any i ∈ [k],
j=1
k X
πij = νj for any j ∈ [kv ],
i=1
e−ε/2 mj ≤ νj ≤ eε/2 mj for any j ∈ [kv ]. where X = [k], V = [kv ], and Cij = d(xi , vj )p . This LP can be reformulated as a minimum cost flow problem similarly to the optimal transport problem, and thus, can be solved in Õ(k 3 ) time [Pele and Werman, 2009]. See Appendix I for details. Thus, exactly solving the above problem is computationally expensive for large k as in the case of optimal transport problems [Cuturi, 2013]. To overcome this issue, we propose to utilize entropic regularization to approximately compute Mm,ε [µ] efficiently. The entropically regularized optimal transport cost [Cuturi, 2013] is defined as follows: X X OTλ (ν, µ) := inf Cij πij + λ πij log πij , π∈Π(ν,µ)
i,j
i,j
where λ > 0 is the regularization parameter. The entropic regularization allows us to compute the optimal transport cost efficiently using the Sinkhorn algorithm [Sinkhorn and Knopp, 1967, Cuturi, 2013], which is an iterative method for solving the regularized optimal transport problem. On the other hand, our focus is computing the projection of a distribution µ onto the LDP polytope Qm,ε with respect to the Wasserstein distance. If we replace the Wasserstein distance with the regularized transport cost, the approximate projection Mλm,ε [µ] is formulated as: Mλm,ε [µ] := arg min OTλ (ν, µ). ν∈Qm,ε
Since x 7→ x1/p is monotone on R+ , this coincides with the exact Wasserstein projection when λ = 0. Note the range of Mλm,ε is still contained in Qm,ε , and thus, Mλm,ε satisfies ε-LDP. For a fixed base measure m, the next proposition quantifies the approximation error of the entropic objective relative to the exact Wasserstein projection. 6
Proposition 4.1 (Entropic approximation error for fixed base measure). Fix m and ε such that Qm,ε ̸= ∅. For µ ∈ P([k]), let ν ∗ ∈ arg min Wp (ν, µ), ν∈Qm,ε
ν λ := Mλm,ε [µ].
Then, we have 0 ≤ Wp (ν λ , µ) − Wp (ν ∗ , µ) ≤ 2λ log(k)
1/p
.
See Appendix K for the proof. By choosing λ sufficiently small, we can make the approximation error arbitrarily small. Naively, the above optimization problem is a bi-level optimization problem, which can be computationally expensive to solve. However, we can show that the above optimization problem can be reformulated as a KL projection onto the intersection of two convex sets on the probability simplex and projection onto each set can be computed efficiently. o n Pk Pkv v Lemma 4.2. Let ∆k×kv := π ∈ Rk×k : i=1 j=1 πij = 1 , and + Rµ := {π ∈ ∆k×kv : π 1 = µ},
Rm,ε := {π ∈ ∆k×kv : π ⊤ 1 ∈ Qm,ε }. P P Assume Qm,ε is nonempty, equivalently e−ε/2 j mj ≤ 1 ≤ eε/2 j mj . Then, we have π∗ = arg
min
π∈Rµ ∩Rm,ε
KL(π∥K),
Mλm,ε [µ] = π∗⊤ 1,
C where K is the Gibbs kernel (depending on λ) defined as Ki,j = Z1λ exp − λi,j , and Zλ is the normalization constant. Furthermore, projection onto Rµ and Rm,ε can be computed as follows: qj∗ ProjKL πij , Rm,ε (π)ij = P i πij
µi ProjKL πij , Rµ (π)ij = P j πij
P ∗ ∗ where sj := = ProjKL Qm,ε (s) = arg minq∈Qm,ε KL(q∥s). Moreover, qj = i πij and q P min{max{eθ sj , e−ε/2 mj }, eε/2 mj }, where θ is a scalar chosen so that j qj∗ = 1. See Appendix J for the proof. The projection q ∗ onto Qm,ε can be computed efficiently using scalar bisection since qj∗ is a monotone function of θ as discussed in Park et al. [2024]. For simplicity, we assume that q ∗ can be computed exactly in the following analysis. The above problem is a special case of Bregman projection onto the intersection of convex sets, and can be solved using the Dykstra’s algorithm [Benamou et al., 2015] based on alternating projections onto Rµ and Rm,ε . In general, Dykstra’s algorithm requires additional correction terms to ensure convergence to the optimal solution and convergence rate cannot be easily characterized. However, utilizing the special structure of the problem, we can prove that a vanilla iterative Bregman projection algorithm without Dykstra’s correction is sufficient to find the optimal solution, which significantly reduces the computational cost. Furthermore, we can show that the algorithm converges linearly to the optimal solution, which is a stronger convergence guarantee than the general case. We provide the details of the algorithm in Algorithm 1. kv To analyze the convergence of the above algorithm, we define the composite mapping T : R++ →
Π K ⊤ (µ/(Kv)) (t) , where Π(s) := ProjKL = T t (v (0) ) for all Qm,ε (s). Note that v K ⊤ (µ/(Kv)) i (xi /yi ) v t ≥ 0. We equip Rk++ with the Hilbert metric defined by dH (x, y) = log max mini (xi /yi ) . v Rk++ by T (v) :=
Then, we have the following convergence guarantee: Theorem 4.3 (Linear convergence of Algorithm 1). Assume that µi > 0 for all i ∈ [k], mj > 0 for all j ∈ [kv ], and Ki,j > 0 for all i ∈ [k], j ∈ [kv ]. (If some µi or mj vanish, one may restrict K to the rows and columns in supp(µ) × supp(m).) Let τ (·) denote the Birkhoff contraction coefficient [Birkhoff, 1957] τ (A) =
dH (Ax, Ay) , dH (x, y) x,y∈Rkv sup
++
7
Algorithm 1 Projection onto LDP polytope with entropic regularization 1: Input: Target distribution µ, base measure m, privacy parameter ε, regularization parameter λ, cost matrix C. 2: Output: Projected distribution Mλ m,ε [µ]. kv 3: Initialize v (0) ∈ R++ . 4: Compute the Gibbs kernel K as [K]i,j = Z1
λ
C
exp − λi,j .
5: for t = 0, 1, 2, . . . until convergence do 6: 7: 8:
Update u(t) = Kvµ(t) . Set s(t) = K ⊤ u(t) . (t) Compute q (t) = ProjKL Qm,ε (s ) via scalar bisection. q (t)
9: Update v (t+1) = s(t) . 10: end for (t) 11: Return Mλ m,ε [µ] = q .
Figure 1: Results for synthetic experiments. Left: Transformed Figure 2: Results for MovieLens distribution with WPM and KPM. Right: Convergence speed of experiments. Our WPM outperAlgorithm 1. forms the baselines. which is strictly less than 1 for any positive matrix A, and c := τ (K ⊤ )τ (K) ∈ (0, 1). Then for any v initial point v (0) ∈ Rk++ , the sequence {v (t) }∞ t=0 generated by Algorithm 1 satisfies dH (v (t) , v ∗ ) ≤ ct dH (v (0) , v ∗ ), where v ∗ is a fixed point of T (unique up to positive scaling). Moreover, defining µ µ , q ∗ := Π(s∗ ), s(t) := K ⊤ , q (t) := Π(s(t) ), s∗ := K ⊤ (t) Kv ∗ Kv we have v ∗ = q ∗ /s∗ (componentwise) and KL(q (t) ∥q ∗ ) ≤ 2ct+1 dH (v (0) , v ∗ ). See Appendix L for the proof. Theorem 4.3 guarantees linear convergence of both the iterates {v (t) } in Hilbert metric and the output distributions {q (t) } in KL divergence. These results ensure the efficiency of Algorithm 1.
5
Numerical Experiments
In this section, we validate our proposed methods through numerical experiments on synthetic data and the MovieLens dataset. We also provide additional experiments on geographical datasets in Appendix M. 5.1
Synthetic Data
To see the difference of privatized distributions between the WPM (ours) and the KPM baseline, we consider a simple synthetic example where X = V = {0, . . . , k − 1} is a discrete set of k elements 8
on a ring, and the cost matrix C is defined as d(i, j) = min{|i − j|, k − |i − j|}. Unless otherwise specified, we set k = 30, ε = 5, and p = 2. We construct an input distribution µ by sampling from a Dirichlet distribution with concentration parameter 0.1 and normalizing the sample to sum to 1. WPM captures the underlying geometry Fig. 1 (Left) shows an example of the projected distributions of µ with WPM and KPM. While KPM produces a uniformly spread distribution, WPM captures the geometric structure and assigns high probability to regions near the true high-probability areas, and low probability to regions that are not. Sinkhorn projection converges rapidly Fig. 1 (Right) shows the convergence of Algorithm 1 in terms of the Wasserstein distance between the input and the output distributions. We see that our proposed algorithm converges rapidly near the exact projection. For large λ, the algorithm converges faster but the approximation error is larger, while for small λ, the algorithm converges slower but the approximation error is smaller, which is consistent with Proposition 4.1. 5.2
MovieLens Dataset
Following the experimental setup of Zamanlooy et al. [2025], we consider a movie recommendation scenario, where the goal is to privately recommend the next movie based on the user’s past movie ratings. While Zamanlooy et al. [2025] consider the problem of genre recommendation, we consider the more practical problem of movie recommendation, where the support (i.e., the set of candidate movies for recommendation) is essentially different from the support of the input distribution (i.e., movies that the user has rated in the past). In such scenarios, f -divergences are not suitable for measuring the distance between the input and output distributions, while Wasserstein distance can still capture the distance between the distributions by utilizing the underlying geometry of the movie space. As a baseline, we consider the f -divergence based method which first samples from the privatized genre distribution and then selects a movie uniformly at random from the selected genre, which is a natural extension of the method proposed by Zamanlooy et al. [2025]. In addition, we also consider the Exponential Mechanism [McSherry and Talwar, 2007] with a utility function defined as the negative expected distance u(µ, j) = −Ei∼µ [d(i, j)]. We use the MovieLens100k dataset [Harper and Konstan, 2015] and divide movies into train and test sets. We construct the input distribution so that it is proportional to the ratings of the movies the user has rated in the train set, and sample a movie from the output distribution of the mechanism over the candidate movies in the test set. As an underlying metric, we compute embeddings of the movies via matrix factorization, capturing user co-consumption patterns, and use cosine distance between the embeddings. In this experiment, we regard the base measure as a hyperparameter and select the best base measure for each method from a set of uniform measures based on the validation performance. We evaluate each method by computing Hit Rate@1, the probability that the output movie is included in the user’s test data with high ratings (≥ 4). The error bars in the figures represent the standard error. WPM works under support mismatch Fig. 2 shows the performance of WPM and KPM in terms of Hit Rate@1. We see that WPM significantly outperforms KPM across all values of ε, and the performance gap becomes larger as ε increases. This is because the WPM can handle the movie distribution directly even under support mismatch by utilizing the underlying geometry of the movie space, while the KPM baseline method suffers from the support mismatch issue and thus, its performance degrades significantly. Exponential Mechanism also performs worse, which may be due to the fact that it only depends on the mean of the input distribution, while WPM can utilize the full distributional information.
6
Conclusion
In this paper, we proposed a novel framework for differentially private sampling from distributions based on Wasserstein geometry. We provide the minimax optimal mechanism in terms of worst-case utility based on Wasserstein projection and we show that the optimal mechanism can be obtained by solving a convex optimization problem. We also proposed a practical algorithm for computing approximate projections using entropic regularization and provided convergence guarantees. Finally, 9
we demonstrated the effectiveness of our proposed method through numerical experiments on synthetic and real-world datasets. Limitations and Future Work Our main contribution lies in the theoretical framework and analysis of optimal private sampling mechanisms based on Wasserstein geometry. Since our framework is quite general and can be applied to various settings, extensive empirical evaluation in specific applications is left for future work. In addition, our focus in this paper is on the local model of differential privacy, and extending our framework to the centralized model of differential privacy is an interesting direction for future research.
10
References Ravindra K. Ahuja, Thomas L. Magnanti, and James B. Orlin. Network flows : theory, algorithms, and applications. Prentice-Hall, Inc., 1993. Hilal Asi, Vitaly Feldman, and Kunal Talwar. Optimal algorithms for mean estimation under local differential privacy. In International Conference on Machine Learning, pages 1046–1056. PMLR, 2022. Jean-David Benamou, Guillaume Carlier, Marco Cuturi, Luca Nenna, and Gabriel Peyré. Iterative bregman projections for regularized transportation problems. SIAM Journal on Scientific Computing, 37(2):A1111–A1138, 2015. Garrett Birkhoff. Extensions of jentzsch’s theorem. Transactions of the American Mathematical Society, 85(1):219–227, 1957. Konstantinos Chatzikokolakis, Catuscia Palamidessi, and Marco Stronati. Constructing elastic distinguishability metrics for location privacy. arXiv preprint arXiv:1503.00756, 2015. Eunjoon Cho, Seth A Myers, and Jure Leskovec. Friendship and mobility: user movement in locationbased social networks. In Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 1082–1090, 2011. Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems, 26, 2013. Apple Differential Privacy Team. Learning with privacy at scale, 2017. Technical report. John C Duchi, Michael I Jordan, and Martin J Wainwright. Local privacy and statistical minimax rates. In 2013 IEEE 54th annual symposium on foundations of computer science, pages 429–438. IEEE, 2013. Cynthia Dwork and Aaron Roth. The algorithmic foundations of differential privacy. Foundations and trends® in theoretical computer science, 9(3-4):211–487, 2014. Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam Smith. Calibrating noise to sensitivity in private data analysis. In Theory of cryptography conference, pages 265–284. Springer, 2006. Úlfar Erlingsson, Vasyl Pihur, and Aleksandra Korolova. Rappor: Randomized aggregatable privacypreserving ordinal response. In Proceedings of the 2014 ACM SIGSAC conference on computer and communications security, pages 1054–1067, 2014. Vitaly Feldman, Audra McMillan, Satchit Sivakumar, and Kunal Talwar. Instance-optimal private density estimation in the wasserstein distance. Advances in Neural Information Processing Systems, 37:90061–90131, 2024. Badih Ghazi, Ravi Kumar, and Pasin Manurangsi. User-level differentially private learning via correlated sampling. Advances in Neural Information Processing Systems, 34:20172–20184, 2021. Hrad Ghoukasian, Bonwoo Lee, and Shahab Asoodeh. Locally optimal private sampling: Beyond the global minimax. arXiv preprint arXiv:2510.09485, 2025. Anming Gu, Sasidhar Kunapuli, Mark Bun, Edward Chien, and Kristjan Greenewald. Differentially private wasserstein barycenters. arXiv preprint arXiv:2510.03021, 2025. F Maxwell Harper and Joseph A Konstan. The movielens datasets: History and context. Acm transactions on interactive intelligent systems (tiis), 5(4):1–19, 2015. Hisham Husain, Borja Balle, Zac Cranko, and Richard Nock. Local differential privacy for sampling. In International Conference on Artificial Intelligence and Statistics, pages 3404–3413. PMLR, 2020. Peter Kairouz, Sewoong Oh, and Pramod Viswanath. Extremal mechanisms for local differential privacy. Advances in neural information processing systems, 27, 2014. 11
Shiva Prasad Kasiviswanathan, Homin K Lee, Kobbi Nissim, Sofya Raskhodnikova, and Adam Smith. What can we learn privately? SIAM Journal on Computing, 40(3):793–826, 2011. Yehuda Koren, Robert Bell, and Chris Volinsky. Matrix factorization techniques for recommender systems. Computer, 42(8):30–37, 2009. Nam Lê Tien, Amaury Habrard, and Marc Sebban. Differentially private optimal transport: Application to domain adaptation. In IJCAI, pages 2852–2858, 2019. Daniel Levy, Ziteng Sun, Kareem Amin, Satyen Kale, Alex Kulesza, Mehryar Mohri, and Ananda Theertha Suresh. Learning with user-level privacy. Advances in Neural Information Processing Systems, 34:12466–12479, 2021. Frank McSherry and Kunal Talwar. Mechanism design via differential privacy. In 48th Annual IEEE Symposium on Foundations of Computer Science (FOCS’07), pages 94–103. IEEE, 2007. Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean. Efficient estimation of word representations in vector space. arXiv preprint arXiv:1301.3781, 2013. Hyun-Young Park, Shahab Asoodeh, and Si-Hyeon Lee. Exactly minimax-optimal locally differentially private sampling. Advances in Neural Information Processing Systems, 37:10274–10319, 2024. Ofir Pele and Michael Werman. Fast and robust earth mover’s distances. In 2009 IEEE 12th international conference on computer vision, pages 460–467. IEEE, 2009. Alain Rakotomamonjy and Ralaivola Liva. Differentially private sliced wasserstein distance. In International Conference on Machine Learning, pages 8810–8820. PMLR, 2021. Francesco Ricci, Lior Rokach, and Bracha Shapira. Introduction to recommender systems handbook. In Recommender systems handbook, pages 1–35. Springer, 2010. David Rodríguez-Vítores, Clément Lalanne, and Jean-Michel Loubes. Learning with differentially private (sliced) wasserstein gradients. arXiv preprint arXiv:2502.01701, 2025. Shai Shalev-Shwartz. Online learning and online convex optimization. Foundations and Trends® in Machine Learning, 4(2):107–194, 2025. Richard Sinkhorn and Paul Knopp. Concerning nonnegative matrices and doubly stochastic matrices. Pacific Journal of Mathematics, 21(2):343–348, 1967. Cédric Villani et al. Optimal transport: old and new, volume 338. Springer, 2009. Vishnu Vinod, Krishna Pillutla, and Abhradeep Guha Thakurta. Invisibleink: High-utility and low-cost text generation with differential privacy. arXiv preprint arXiv:2507.02974, 2025. Behnoosh Zamanlooy, Mario Diaz, and Shahab Asoodeh. Locally private sampling with public data. In International Conference on Artificial Intelligence and Statistics, pages 622–630. PMLR, 2025.
A
Detailed Discussion of Related Work
In this section, we show that our framework can be viewed as a generalization of the classical estimation problem in DP. Let us consider the matched-support setting X = V where the user has a single data point x ∈ X , and the goal is to design a randomized mechanism L : X → X that outputs a privatized version of x while satisfying ε-LDP. The utility of the mechanism is typically measured by the mean squared error E[∥x − L(x)∥22 ]. On the other hand, in our framework, let us consider the case where the input distribution µ is a Dirac measure at x, i.e., µ = δx . If we set d(x, y) = ∥x − y∥2 and p = 2, then the Wasserstein distance between M[δx ] and δx is given by Z W22 (M[δx ], δx ) = ∥x − y∥22 dM[δx ](y) = E[∥x − L(x)∥22 ]. X
12
This matches the mean squared error of the mechanism L in the classical estimation problem. Therefore, classical estimation under LDP can be seen as a special case of our framework where the input distribution is a Dirac measure. Note that estimation under LDP typically focuses on estimating a mean across users, and thus it poses an unbiasedness constraint on the mechanism. In such a case, the problem is different from our framework, which does not require unbiasedness.
B
Auxiliary Results
Lemma B.1. A mechanism M satisfies ε-LDP if and only if there exists a base measure m such that for all µ ∈ P(X ), M[µ] ∈ Qm,ε . Proof. If Part. Suppose there exists a base measure m such that for all µ ∈ P(X ), M[µ] ∈ Qm,ε . Then, we have for all µ, µ′ ∈ P(X ) and for all measurable subsets A ⊆ V: M[µ](A) ≤ eε/2 m(A) M[µ′ ](A) ≥ e−ε/2 m(A) and thus, M[µ](A) ≤ eε M[µ′ ](A). This shows that M satisfies ε-LDP. Only if Part. Suppose M satisfies ε-LDP. Fix an arbitrary µ0 ∈ P(X ) and define a reference measure λ = M[µ0 ]. For any µ ∈ P(X ) and measurable A ⊆ V, ε-LDP gives M[µ](A) ≤ eε λ(A). Hence M[µ] ≪ λ, so by the Radon-Nikodym theorem fµ := dM[µ] exists for each µ. dλ Moreover, for any µ, µ′ ∈ P(X ), ε-LDP implies M[µ](A) ≤ eε M[µ′ ](A) for all measurable A. Taking its Radon-Nikodym derivative with respect to λ yields fµ (x) ≤ eε fµ′ (x)
for λ-almost every x.
(2)
Define f := ess supµ∈P(X ) fµ ,
f := ess inf µ∈P(X ) fµ .
Since λ is a probability measure, especially a σ-finite measure, there is a countable set C ⊂ P(X ) such that f (x) = sup fµ (x),
f (x) = inf fµ (x)
µ∈C
µ∈C
for λ-almost every x. From Eq. (2), we have f (x) ≤ eε f (x) for λ-almost every x. We q now construct the base measure m by defining its density with respect to λ as dm(x) = f (x)f (x) dλ(x). For any µ ∈ P(X ), we have f (x) ≤ fµ (x) ≤ f (x) λ-almost everywhere. Thus, λ-almost everywhere, q q q fµ (x) ≤ f (x) · f (x) ≤ eε f (x) · f (x) = eε/2 f (x)f (x) q q q fµ (x) ≥ f (x) · f (x) ≥ e−ε f (x) · f (x) = e−ε/2 f (x)f (x). Integrating these bounds over any measurable subset A ⊆ V with respect to λ, we obtain: e−ε/2 m(A) ≤ M[µ](A) ≤ eε/2 m(A). Therefore, M[µ] ∈ Qm,ε for all µ ∈ P(X ). 13
Proposition B.2. For x ∈ X , let δx be the Dirac measure at x. Then, if M is a Wasserstein projection, we have sup Wp (M[µ], µ) = sup Wp (M[δx ], δx ). x∈X
µ∈P(X )
Proof. Let Q := Qm,ε and assume that M[µ] ∈ arg min Wp (ν, µ) ν∈Q
for every µ ∈ P(X ). The inequality sup Wp (M[δx ], δx ) ≤ x∈X
sup Wp (M[µ], µ) µ∈P(X )
is immediate since δx ∈ P(X ). For the reverse inequality, fix µ ∈ P(X ) and define Z νx := M[δx ],
ν̄ :=
νx dµ(x). X
Since Q is convex, ν̄ ∈ Q. By the projection property, Wpp (M[µ], µ) ≤ Wpp (ν̄, µ). Let πx (dy, dz) = νx (dy)δx (dz) and define Z π :=
πx dµ(x). X
Then π is a coupling of ν̄ and µ. Hence Wpp (ν̄, µ) ≤
Z
d(y, z)p dπ(y, z)
X ×X
Z Z = ZX =
d(y, x)p dνx (y) dµ(x)
V
Wpp (νx , δx ) dµ(x)
X
Z =
Wpp (M[δx ], δx ) dµ(x)
X
≤ sup Wpp (M[δx ], δx ). x∈X
Therefore Wp (M[µ], µ) ≤ sup Wp (M[δx ], δx ). x∈X
Taking supremum over µ completes the proof.
C
Proof of Proposition 3.2
Proof. Fix µ, µ′ ∈ P(X ) and A ⊆ V measurable. By the definition of Qm,ε , Mm,ε [µ](A) ≤ βm(A) Mm,ε [µ′ ](A) ≥ αm(A) where α := e−ε/2 and β := eε/2 . Hence Mm,ε [µ](A) ≤ eε Mm,ε [µ′ ](A). This proves that Mm,ε satisfies ε-LDP. 14
D
Proof of Proposition 3.3
Proof. For any ε-LDP mechanism M, by Lemma B.1, there exists a base measure m such that for all µ, M[µ] ∈ Qm,ε . By the definition of Mm,ε , Mm,ε [µ] is the projection of µ onto Qm,ε with respect to Wp . From the property of projection, for all ν ∈ Qm,ε , Wp (Mm,ε [µ], µ) ≤ Wp (ν, µ). Taking ν = M[µ] ∈ Qm,ε gives the desired result.
E
Proof of Theorem 3.4
Proof. By Theorem 3.3, for any ε-LDP mechanism M, there exists a base measure m such that for all µ, Wp (M[µ], µ) ≥ Wp (Mm,ε [µ], µ). Hence inf
sup Wp (M[µ], µ) ≥ inf sup Wp (Mm,ε [µ], µ) m
ε-LDP M µ
µ
≥ inf max ϕi (m)1/p . m
i
The reverse inequality is immediate since Mm,ε is ε-LDP for every m. Thus, it suffices to show that the infimum over m is attained by some m∗ . By Step 4 in the proof of Theorem 3.5, m 7→ maxi ϕi (m) is a convex function on Rk+v . Moreover, the P constraint set {m ∈ Rk+v : 1/β ≤ j mj ≤ 1/α} is a compact convex set. Therefore, the infimum is attained by some m∗ .
F
Proof of Proposition 3.5
Proof. Throughout, let X = [k], V = [kv ], α := e−ε/2 , and β := eε/2 . Step 1: Coordinate form of Qm,ε . In the discrete case, for any nonnegative vector m ∈ Rk+v and any probability vector ν ∈ P([kv ]), ν ∈ Qm,ε ⇐⇒ αmi ≤ νi ≤ βmi
(∀i ∈ [kv ]).
The “only if” direction follows by taking A = {i} in the definition of Qm,ε . The converse follows by summing coordinate-wise inequalities over any A ⊆ [kv ]. Step 2: Reduction to Dirac inputs.
By Proposition B.2, for projection mechanisms, sup Wp Mm,ε [µ], µ = max Wp Mm,ε [δi ], δi . i∈[k]
µ∈P(X )
Since x 7→ xp is strictly increasing on R+ , minimizing the left-hand side over m is equivalent to minimizing max Wpp Mm,ε [δi ], δi . i∈[k]
Step 3: Inner optimization for each i. Fix i ∈ [k]. For any ν ∈ P([kv ]), since one marginal is δi , the only coupling of (ν, δi ) is π(j, ℓ) = νj 1{ℓ=i} . 15
Therefore, Wpp (ν, δi ) =
kv X
d(vj , xi )p νj =
j=1
kv X
Cij νj .
j=1
Hence Wpp Mm,ε [δi ], δi = min
kv X
q∈Rkv
Cij qj :
j=1
kv X
qj = 1, qj ≥ 0, αmj ≤ qj ≤ βmj (∀j)
.
j=1
The above problem is a fractional knapsack problem: for each i ∈ [k], sort indices by increasing Cij and greedily fill capacities (β − α)mj . Therefore, ϕi (m) can be expressed in the following closed form: X X ϕi (m) = α Cij mj + β Cij mj j:Cij >τi (m)
j:Cij <τi (m)
X
+ τi (m) 1 − α
mj − β
j:Cij >τi (m)
= τi (m) +
kv X
X
mj
j:Cij <τi (m)
(α(Cij − τi (m))+ − β(τi (m) − Cij )+ )mj ,
j=1
where τi (m) is defined as τi (m) := min{t :
Pkv
j=1 (βmj 1{Cij ≤t} + αmj 1{Cij >t} ) ≥ 1}.
Combining Steps 2 and 3, !1/p inf
sup Wp Mm,ε [µ], µ =
v m∈Rk + µ∈P(X )
inf max ϕi (m)
v i∈[k] m∈Rk +
.
Therefore, minimizing worst-case utility is equivalent to min max ϕi (m).
v i∈[k] m∈Rk +
Step 4: Convexity in the base measure. Define kv X qj = 1, αmj ≤ qj ≤ βmj (∀j) . Q(m) := q ∈ Rk+v : j=1
For m, m′ ∈ Rk+v and θ ∈ [0, 1], let mθ := θm + (1 − θ)m′ . If either Q(m) or Q(m′ ) is empty, convexity of the extended-value function ϕi is immediate. Otherwise, for any q ∈ Q(m) and q ′ ∈ Q(m′ ), qθ := θq + (1 − θ)q ′ satisfies X qθ,j = 1, αmθ,j ≤ qθ,j ≤ βmθ,j (∀j), j
so qθ ∈ Q(mθ ). Hence ϕi (mθ ) ≤
kv X j=1
Cij qθ,j = θ
kv X j=1
Cij qj + (1 − θ)
kv X
Cij qj′ .
j=1
Taking infimum over q ∈ Q(m) and q ′ ∈ Q(m′ ) gives ϕi θm + (1 − θ)m′ ≤ θϕi (m) + (1 − θ)ϕi (m′ ); Thus each ϕi is convex (as an extended-value function on Rk+v ), and therefore m 7→ maxi ϕi (m) is convex. This completes the proof. 16
G
Algorithm details for base-measure optimization
Define f (m) := max ϕi (m),
Dp :=
i∈[k]
max i∈[k], j∈[kv ]
Cij .
Let kv X 1 1 . M := m ∈ Rk+v : ≤ mj ≤ β α j=1
By Step 3 in the proof of Theorem 3.5, ϕi (m) has an equivalent LP form. For fixed i and m, write qj = αmj + rj . Then kv kv kv kv X X X X ϕi (m) = α Cij mj + min Cij rj : 0 ≤ rj ≤ (β − α)mj (∀j), rj = 1 − α mj . r j=1
j=1
j=1
j=1
Hence feasibility is equivalent to 0≤1−α
X
mj ≤ (β − α)
j
X
mj ,
j
i.e., kv X 1 1 mj ≤ . ≤ β α j=1
Let τi (m) denote the threshold cost of the greedy solution. To obtain a subgradient with respect to m, we derive the LP dual. Returning to the original primal form (Step 3), introduce dual variables η ∈ R P for j qj = 1, sj ≥ 0 for qj ≥ αmj , and tj ≥ 0 for qj ≤ βmj . The Lagrangian is X P P P L= Cij qj + η 1 − j qj − j sj (qj − αmj ) − j tj (βmj − qj ) j
=η+
X X (αsj − βtj )mj + (Cij − η − sj + tj ) qj . | {z } j j =0 at optimum
Taking inf q is finite iff sj − tj = Cij − η for all j. Hence the dual is ϕi (m) =
max
v η∈R, s,t∈Rk +
kv o n X η+ (αsj − βtj )mj : sj − tj = Cij − η (∀j) .
∗
j=1 ∗
∗
Hence, for any dual optimizer (η , s , t ), Danskin’s theorem gives (i)
gj (m) := αs∗j − βt∗j ,
j ∈ [kv ],
and g (i) (m) ∈ ∂ϕi (m). By complementary slackness, t∗j > 0 (i.e. η ∗ > Cij ) implies qj = βmj , while s∗j > 0 (i.e. η ∗ < Cij ) implies qj = αmj . Hence η ∗ coincides with the threshold cost τi (m) of the greedy solution in Step 3. Setting τi (m) := η ∗ , the constraint sj − tj = Cij − τi (m) with sj , tj ≥ 0 is uniquely decomposed as s∗j = (Cij − τi (m))+ ,
t∗j = (τi (m) − Cij )+ ,
and we obtain the explicit subgradient (i)
gj (m) = α(Cij − τi (m))+ − β(τi (m) − Cij )+ , Therefore (i)
|gj (m)| ≤ βDp ,
∥g (i) (m)∥∞ ≤ βDp . 17
j ∈ [kv ].
∗
If i∗ ∈ arg maxi ϕi (m), then g (i ) (m) ∈ ∂f (m). Consequently, for any m, m′ ∈ M, |ϕi (m) − ϕi (m′ )| ≤ βDp ∥m − m′ ∥1 ,
|f (m) − f (m′ )| ≤ βDp ∥m − m′ ∥1 .
We optimize f directly over M using generalized-KL mirror descent. At iteration t, choose it ∈ arg maxi ϕi (m(t) ) and let g (t) ∈ ∂ϕit (m(t) ) ⊆ ∂f (m(t) ). With stepsize ηt > 0, perform (t+1)
m̃i
(t)
(t)
= mi exp −ηt gi
,
i ∈ [kv ],
then project by generalized KL: m(t+1) = arg min DgKL m m̃(t+1) , m∈M
where DgKL (u∥v) :=
kv X i=1
ui ui log − ui + v i . vi
P (t+1) P This projection has a closed form. Let St := i m̃i . Then for fixed total mass s = i mi , KKT gives the minimizer m = (s/St )m̃(t+1) , and minimizing over s ∈ [1/β, 1/α] reduces to clipping St to this interval. (1/β) (t+1) m̃ , St < 1/β, St m(t+1) = m̃(t+1) , 1/β ≤ St ≤ 1/α, (1/α) (t+1) m̃ , St > 1/α. St Return T
1 X (t) m . T t=1
m̄T :=
Proposition G.1. Let m∗ ∈ arg minm∈M f (m) and initialize m(1) =
1 1. αkv
Choose a constant stepsize p ηt ≡ η :=
2(1 + log k) √ . βDp T
Then βDp f (m̄T ) − f (m ) ≤ α
r
∗
2(1 + log k) β =O Dp T α
r
log k T
! .
Proof. By construction, g (t) ∈ ∂f (m(t) ) for every t. By the regret bound for the unnormalized exponentiated gradient algorithm [Shalev-Shwartz, 2025, Theorem 2.23], for any m ∈ M, T T kv X DgKL (m∥m(1) ) η X X (t) (t) 2 . ⟨g (t) , m(t) − m⟩ ≤ + mi gi η 2 t=1 t=1 i=1
Using ∥g (t) ∥∞ ≤ βDp and
(t) i mi ≤ 1/α,
P
kv X
(t)
mi
(t) 2
gi
i=1
18
≤
(βDp )2 . α
Next, fix any m ∈ M. Since m(1) = (1/(αkv ))1, DgKL (m∥m(1) ) =
kv X
mj log(αkv mj ) −
j=1
Let pj := αmj , so
kv X
mj +
j=1
1 . α
P
j pj ≤ 1. Then k
kv X
v log k 1X pj log(kpj ) ≤ , mj log(αkv mj ) = α α j=1 j=1
because
P
j pj log pj ≤ 0 and
P
j pj ≤ 1. Therefore
DgKL (m∥m(1) ) ≤
1 + log k . α
Applying the regret bound with m = m∗ and dividing by T gives T
1 X (t) (t) 1 + log k η(βDp )2 ⟨g , m − m∗ ⟩ ≤ + . T t=1 αηT 2α By convexity of f , f (m(t) ) − f (m∗ ) ≤ ⟨g (t) , m(t) − m∗ ⟩. PT Using Jensen’s inequality for m̄T = T1 t=1 m(t) , f (m̄T ) − f (m∗ ) ≤
T 1 + log k η(βDp )2 1X f (m(t) ) − f (m∗ ) ≤ + . T t=1 αηT 2α
Substituting p η=
2(1 + log k) √ βDp T
yields βDp f (m̄T ) − f (m ) ≤ α ∗
r
2(1 + log k) . T
v For computation, first sort {Cij }kj=1 for each i ∈ [k] once in O(k 2 log k) time. Then each iteration evaluates all ϕi and one active subgradient in O(k 2 ) time, so with ! (β/α)2 Dp2 log k T =O ε2opt
iterations, the total complexity is (β/α)2 Dp2 log k O k log k + k ε2opt 2
H
2
! .
Optimal Base Measure on the Sphere
As an example of a continuous setting, we consider X = S d := {x ∈ Rd+1 : ∥x∥2 = 1} equipped with the Euclidean distance d(x, y) = ∥x − y∥2 . Utilizing the rotational symmetry of the sphere, we can derive the optimal base measure analytically. 19
Proposition H.1 (Optimal base measure on the sphere). Fix p ≥ 1, and let σ denote the uniform probability measure on S d . Then, the optimal base measure is given by m∗ = α∗ σ, where α∗ is a scale parameter defined as α∗ =
1 , e−ε/2 + eε/2 − e−ε/2 Sd (t∗ )
with d−2
Γ((d + 1)/2) Cd = √ , πΓ(d/2) Z 1 Sd (t) = fd (u) du,
fd (u) = Cd (1 − u2 ) 2 , gp (u) = (2 − 2u)p/2 ,
t
Z 1 Ud,p (t) =
Z 1 gp (u)fd (u) du,
Hd,p =
gp (u)fd (u) du, −1
t
and t∗ ∈ (−1, 1) is the unique solution of eε/2 − e−ε/2 Ud,p (t∗ ) + e−ε/2 Hd,p = gp (t∗ ) e−ε/2 + eε/2 − e−ε/2 Sd (t∗ ) . Moreover, the map p 7→ α∗ (p) is decreasing. The optimal mechanism Mm∗ ,ε can be interpreted as a generalization of PrivUnit [Asi et al., 2022], which is an optimal LDP mechanism for privatizing unit vectors on the sphere, except that we do not require the unbiasedness condition since our focus is not on mean estimation. Interestingly, the scale parameter α∗ is decreasing in p: as p increases, the cost of transporting mass over long distances becomes more expensive, and thus it is preferable to reduce the minimum amount of mass (αe−ε/2 σ) that must be assigned to locations far away. Proof. We divide the proof into four steps. Step 1: Reduction to rotation-invariant base measures. Let σ denote the uniform probability measure on S d . By Proposition B.2, it is enough to minimize F (m) := sup Wp δx , ΠQm,ε (δx ) . x∈S d
For u ∈ S d , define ϕu (m) :=
inf
ν∈Qm,ε
Wpp (δu , ν).
Since the source is a Dirac mass, Wpp (δu , ν) =
Z Sd
∥u − y∥p2 dν(y),
so the optimization defining ϕu (m) is linear over Qm,ε . Fix, for each u ∈ S d , an optimizer νu ∈ Qm,ε , and define Z Z m̄(A) := (R# m)(A) dλ(R), ν̄x (A) := (R# νR−1 x )(A) dλ(R). SO(d+1)
SO(d+1) d
For every measurable A ⊂ S , e−ε/2 m̄(A) ≤ ν̄x (A) ≤ eε/2 m̄(A), hence ν̄x ∈ Qm̄,ε . Therefore, Z ϕx (m̄) ≤ Sd
∥x − y∥p2 dν̄x (y).
20
Using rotational invariance of the Euclidean norm, Z Z Z ∥x − y∥p2 dν̄x (y) = ∥x − y∥p2 d(R# νR−1 x )(y) dλ(R) Sd SO(d+1) S d Z Z = ∥R−1 x − z∥p2 dνR−1 x (z) dλ(R) d SO(d+1) S Z = ϕR−1 x (m) dλ(R) SO(d+1)
≤ sup ϕu (m). u∈S d
Hence sup ϕx (m̄) ≤ sup ϕu (m). u∈S d
x∈S d
Since F (m)p = supx∈S d ϕx (m), we obtain F (m̄) ≤ F (m). Thus an optimal base measure can be taken rotation-invariant, i.e., m = ασ for some α > 0. Step 2: Structure of the optimal projection. Fix x ∈ S d and α > 0. Any ν ∈ Qασ,ε is absolutely continuous with respect to σ with density q = dν/dσ satisfying αe−ε/2 ≤ q(y) ≤ αeε/2 for σ-almost every y, and Z q(y) dσ(y) = 1. Sd
Since Wpp (δx , ν) =
Z Sd
∥x − y∥p2 q(y) dσ(y),
minimizing the Wasserstein cost is equivalent to minimizing the above linear functional. Let u = ⟨x, y⟩ ∈ [−1, 1] and define gp (u) := (2 − 2u)p/2 . Then ∥x − y∥p2 = gp (⟨x, y⟩), and gp is strictly decreasing on [−1, 1]. By the bathtub principle, the minimizer saturates the upper bound on a spherical cap and the lower bound on its complement. Thus, for some threshold t ∈ [−1, 1], qx,t (y) = αeε/2 1{⟨x,y⟩≥t} + αe−ε/2 1{⟨x,y⟩<t} . Let Sd (t) := σ {y ∈ S d : ⟨x, y⟩ ≥ t} . The normalization condition
R
qx,t dσ = 1 gives αeε/2 Sd (t) + αe−ε/2 (1 − Sd (t)) = 1,
hence α=
1 e−ε/2 +
eε/2 − e−ε/2 21
Sd (t)
.
Step 3: Reduction to a one-dimensional objective. Under the uniform measure σ, the random variable u = ⟨x, y⟩ has density Γ((d + 1)/2) Cd = √ . πΓ(d/2)
d−2
fd (u) = Cd (1 − u2 ) 2 , Hence Z 1
fd (u) du.
Sd (t) = t
Define Z 1 Ud,p (t) :=
gp (u)fd (u) du, t
Z 1 Hd,p :=
gp (u)fd (u) du. −1
For νx,t induced by qx,t , Wpp (δx , νx,t ) = αeε/2 Ud,p (t) + αe−ε/2 Hd,p − Ud,p (t) = α eε/2 − e−ε/2 Ud,p (t) + e−ε/2 Hd,p . Let A := eε/2 − e−ε/2 ,
B := e−ε/2 .
Substituting α = 1/(B + ASd (t)), we get the one-dimensional objective J(t) := Wpp (δx , νx,t ) =
AUd,p (t) + BHd,p . B + ASd (t)
Step 4: Optimality condition and uniqueness. Write U (t) := Ud,p (t), S(t) := Sd (t), H := Hd,p , and g(t) := gp (t). Since U ′ (t) = −g(t)fd (t),
S ′ (t) = −fd (t),
differentiation yields J ′ (t) =
Afd (t) G(t), (B + AS(t))2
where G(t) := AU (t) + BH − g(t) B + AS(t) . Also, G′ (t) = AU ′ (t) − g ′ (t) B + AS(t) − g(t)AS ′ (t) = −g ′ (t) B + AS(t) . Because g is strictly decreasing, g ′ (t) < 0 on (−1, 1); and B + AS(t) > 0. Therefore G′ (t) > 0 on (−1, 1), so G is strictly increasing. At t = −1 and t = 1, G(−1) = (A + B)H − g(−1)(A + B) = (A + B) H − 2p < 0, G(1) = BH > 0, p
since g(−1) = 2 and H < 2p . Thus, by continuity and strict monotonicity, there exists a unique t∗ ∈ (−1, 1) such that G(t∗ ) = 0, i.e. eε/2 − e−ε/2 Ud,p (t∗ ) + e−ε/2 Hd,p = gp (t∗ ) e−ε/2 + eε/2 − e−ε/2 Sd (t∗ ) .
22
Since A > 0, fd (t) > 0 on (−1, 1), and (B + AS(t))2 > 0, we have sign(J ′ (t)) = sign(G(t)). Hence J is strictly decreasing on (−1, t∗ ) and strictly increasing on (t∗ , 1). Therefore t∗ is the unique global minimizer of J, equivalently of Wpp (δx , νx,t ). Finally, substituting t∗ into the normalization formula gives α∗ =
1 . e−ε/2 + eε/2 − e−ε/2 Sd (t∗ )
Thus m∗ = α∗ σ is optimal, with t∗ uniquely characterized by the stated scalar equation. Step 5: Monotonicity of p 7→ α∗ (p). From the formula for α∗ and the monotonicity of Sd (t), it suffices to show that p 7→ t∗ (p) is decreasing on [1, ∞). Write √ c± := e±ε/2 , r(u) := 2 − 2u, so that gp (u) = r(u)p . For t ∈ (−1, 1), define Z(t) := c− + (c+ − c− )Sd (t), and the probability measure c+ 1{u≥t} + c− 1{u<t} fd (u) du µt (du) := . Z(t) For each p ≥ 1, let t∗ (p) be the unique solution from Step 4. Its defining equation is equivalent to Z 1 gp (u) µt∗ (p) (du) = gp (t∗ (p)), −1
equivalently, 1/p = r t∗ (p) . r(u) µt∗ (p) (du)
Z 1
p
−1
Define Z 1 Mp (t) :=
1/p , r(u)p µt (du)
Fp (t) := Mp (t) − r(t).
−1
Then t∗ (p) is the unique zero of Fp . Fix t ∈ (−1, 1). Since c± > 0, fd (u) > 0 for u ∈ (−1, 1), and Z(t) > 0, µt has positive density on (−1, 1). As r is positive and nonconstant on (−1, 1), strict monotonicity of power means implies p1 < p2
=⇒
Mp1 (t) < Mp2 (t).
Let ti := t∗ (pi ) for 1 ≤ p1 < p2 . Since Mp1 (t1 ) = r(t1 ), Mp2 (t1 ) > Mp1 (t1 ) = r(t1 ), hence Fp2 (t1 ) > 0. Also, Fp is continuous in t, and by dominated convergence, 1/p
lim Fp (t) = Hd,p − 2 < 0,
t↓−1
1/p
lim Fp (t) = Hd,p > 0. t↑1
1/p
(Indeed, r(u) ≤ 2 and r(u) < 2 for u > −1, so Hd,p < 2.) Since Fp2 has a unique zero, Fp2 (t1 ) > 0 implies t∗ (p2 ) < t∗ (p1 ). Thus p 7→ t∗ (p) is strictly decreasing on [1, ∞). This completes the proof. 23
I
LP Reformulation as Minimum-Cost Flow
We show that the LP formulation of the Wasserstein projection in Section 4 can be reformulated as a minimum-cost flow problem on a bipartite network. Recall the LP: kv k X X
min
k×kv v ν∈Rk + , π∈R+ i=1 j=1
subject to
kv X
Cij πij
πij = µi for any i ∈ [k],
j=1
k X
πij = νj for any j ∈ [kv ],
i=1
e−ε/2 mj ≤ νj ≤ eε/2 mj for any j ∈ [kv ]. We construct a directed graph G = (V, E) as follows. Nodes. The node set is V = {s} ∪ {ui : i ∈ [k]} ∪ {vj : j ∈ [kv ]} ∪ {t}, where s is the source and t is the sink. Edges and capacities. • Supply edges: For each i ∈ [k], add edge (s, ui ) with capacity [µi , µi ] (i.e., the flow is fixed to µi ) and cost 0. • Transport edges: For each i ∈ [k] and j ∈ [kv ], add edge (ui , vj ) with capacity [0, µi ] and cost Cij . • Demand edges: For each j ∈ [kv ], add edge (vj , t) with capacity [e−ε/2 mj , eε/2 mj ] and cost 0. P The total supply at s equals i µi = 1, and the demand edges enforce that the total flow into t also P P equals 1 (since Qm,ε ̸= ∅ implies e−ε/2 j mj ≤ 1 ≤ eε/2 j mj ). Correspondence. Let f be a feasible flow on G. Setting πijP= f (ui , vj ) and νj = f (vj , t) recovers a feasible solution of the LP, and vice versa. The objective i,j Cij f (ui , vj ) coincides with the LP objective. Thus the LP is equivalent to a minimum-cost flow problem on G. Complexity. The network has |V | = k + kv + 2 nodes and |E| = k + k · kv + kv edges. This is analogous to the standard network-flow formulation of the discrete optimal transport problem, where the only difference is that the demand edges carry box constraints from the LDP polytope instead of fixed demands. Via the successive shortest path algorithm, the minimum-cost flow problem can be solved in Õ(k 3 ) time [Ahuja et al., 1993].
J
Proof of Lemma 4.2
v Proof. Let C ∈ Rk×k be the cost matrix. Recall the entropically regularized transport cost: + X X OTλ (µ, ν) = min Cij πij + λ πij log πij .
π∈∆k×kv i,j π 1=µ, π ⊤ 1=ν
i,j
Hence Mλm,ε [µ] = arg min OTλ (µ, ν) ν∈Qm,ε X X = arg min Cij πij + λ πij log πij . π∈∆k×kv i,j π 1=µ, π ⊤ 1∈Qm,ε
24
i,j
Since log Kij = −
Cij − log Zλ , λ
we have KL(π∥K) =
X
=
X
πij log
i,j
πij Kij
πij log πij +
i,j
X 1X Cij πij + (log Zλ ) πij . λ i,j i,j
On the constraint set π 1 = µ, the total mass is fixed and equals X X πij = µi = 1. i,j
i
Cij πij + λ
X
Therefore minimizing X i,j
πij log πij
i,j
is equivalent (up to constants independent of π) to minimizing KL(π∥K). This proves π∗ = arg
min
π∈Rµ ∩Rm,ε
KL(π∥K).
By construction, ν∗ := π∗ ⊤ 1 belongs to Qm,ε and attains min OTλ (µ, ν),
ν∈Qm,ε
hence Mλm,ε [µ] = ν∗ . Next we derive the closed-form KL projections. Projection onto Rµ : PRµ (π) = arg min KL(π̂∥π)
X
s.t.
π̂≥0
π̂ij = µi (∀i).
j
The Lagrangian optimality condition gives π̂ij log + ai = 0, πij so π̂ij = πij e−ai . Enforcing row sums, µi , e−ai = P j πij and therefore µi PRµ (π)ij = P πij . j πij Projection onto Rm,ε :
write sj :=
P
i πij . The constraint π
X
e−ε/2 mj ≤ qj ≤ eε/2 mj ,
⊤
1 ∈ Qm,ε is
qj = 1,
(∀j),
j
where qj :=
P
i π̂ij . For fixed q, we have the decomposition
KL(π̂∥π) =
X
π̂ij πij
π̂ij qj πij /sj
π̂ij log
i,j
=
X i,j
π̂ij log
+
X j
25
qj log
qj sj
≥
X j
qj log
qj sj
,
where equality holds iff π̂ij = qj πij /sj . Therefore, the column-marginal projection reduces to X qj ∗ . q = arg min qj log q∈Qm,ε sj j Let α := e−ε/2 and β := eε/2 . The KKT conditions for X qj min qj log q sj j X s.t. qj = 1, αmj ≤ qj ≤ βmj (∀j) j
yield qj∗ = min{max{eθ sj , αmj }, βmj }, P ∗ where θ ∈ R is chosen so that j qj = 1. Define X ψ(θ) := min{max{eθ sj , αmj }, βmj }. j
Then ψ is continuous and nondecreasing, with X lim ψ(θ) = α mj , θ→−∞
Since Qm,ε is nonempty, α ψ(θ) = 1. Finally,
lim ψ(θ) = β
θ→+∞
j
P
j mj
≤ 1 ≤ β
X
mj .
j
P
j mj , hence there exists at least one θ such that
qj∗ PRm,ε (π)ij = P πij . i πij This completes the proof.
K
Proof of Proposition 4.1
Proof. Let Q := Qm,ε and N := kkv . P For any ν ∈ P([kv ]) and any coupling π ∈ Π(ν, µ) with i,j πij = 1, we have X − log N ≤ πij log πij ≤ 0, i,j
where 0 log 0 := 0. Therefore, X X X X Cij πij − λ log N ≤ Cij πij + λ πij log πij ≤ Cij πij . i,j
i,j
i,j
i,j
Taking the minimum over π ∈ Π(ν, µ) yields, for every ν, OT(ν, µ) − λ log N ≤ OTλ (ν, µ) ≤ OT(ν, µ). Let ν ∗ ∈ arg min Wp (ν, µ), ν∈Q
ν λ := Mλm,ε [µ].
Since x 7→ x1/p is strictly increasing on R+ , we also have ν ∗ ∈ arg min OT(ν, µ). ν∈Q
By definition of Mλm,ε [µ], ν λ ∈ arg min OTλ (ν, µ). ν∈Q
26
(3)
Using (3), OT(ν λ , µ) ≤ OTλ (ν λ , µ) + λ log N ≤ OTλ (ν ∗ , µ) + λ log N ≤ OT(ν ∗ , µ) + λ log N. On the other hand, by optimality of ν ∗ for OT over Q, OT(ν ∗ , µ) ≤ OT(ν λ , µ). Hence 0 ≤ OT(ν λ , µ) − OT(ν ∗ , µ) ≤ λ log N. Finally, let b := Wpp (ν ∗ , µ),
a := Wpp (ν λ , µ),
so 0 ≤ a − b ≤ λ log N . Because a = b + (a − b) and x 7→ x1/p is subadditive on R+ for p ≥ 1, Wp (ν λ , µ) − Wp (ν ∗ , µ) = a1/p − b1/p ≤ (a − b)1/p ≤ (λ log N )1/p . This proves the claim.
L
Proof of Theorem 4.3
First, we show some technical lemmas. Lemma L.1 (Nonexpansiveness of simplex-projection ratio). Assume mj > 0 for all j ∈ [kv ], and define Π(x) := arg min KL(q∥x), q∈Qm,ε
F (x) :=
Π(x) x
v v , . Then for all x, y ∈ Rk++ for x ∈ Rk++
dH (F (x), F (y)) ≤ dH (x, y). Proof. Let ℓi := e−ε/2 mi and ui := eε/2 mi . From the definition of Π, we have log F (x)i = log Π(x)i − log xi = clip(log xi + θ(x), [log ℓi , log ui ]) − log xi = clip(log xi + θ(x), [log ℓi , log ui ]) − log xi − θ(x) + θ(x) | {z } =:a(x)i
= ϕi (a(x)i ) − a(x)i + θ(x) = ψi (a(x)i ) + θ(x), where ϕi (z) := clip(z, [log ℓi , log ui ]) is monotonically nondecreasing and 1-Lipschitz, and ψi (z) := ϕi (z) − z is monotonically nonincreasing and 1-Lipschitz. Let m = mini {a(x)i − a(y)i } and M = maxi {a(x)i − a(y)i }. From the invariance of Hilbert’s metric under scaling, we have dH (x, y) = dH (exp(log x) · exp(θ(x)), exp(log y) · exp(θ(y))) = dH (exp(a(x)), exp(a(y))) = M − m, dH (F (x), F (y)) = dH (exp(ψ(a(x)) + θ(x)), exp(ψ(a(y)) + θ(y))) = dH (exp(ψ(a(x))), exp(ψ(a(y)))). Thus, it suffices to show that ψi (a(x)i ) − ψi (a(y)i ) ∈ [−M, −m]. 27
Case 1: m ≤ 0 ≤ M . If a(x)i ≥ a(y)i , we have ψi (a(x)i ) − ψi (a(y)i ) ≤ 0 ≤ −m since ψi is nonincreasing. In addition, since ψi is 1-Lipschitz, ψi (a(x)i ) − ψi (a(y)i ) ≥ −|a(x)i − a(y)i | = −(a(x)i − a(y)i ) ≥ −M . If a(x)i < a(y)i , we have ψi (a(x)i ) − ψi (a(y)i ) ≥ 0 since ψi is nonincreasing, and ψi (a(x)i ) − ψi (a(y)i ) ≤ |a(x)i − a(y)i | = a(y)i − a(x)i ≤ −m. Case 2: 0 < m ≤ M . In this case, a(x)i is always P greater than P a(y)i . From the monotonicity of ϕi , we have ϕ(a(x))i ≥ ϕ(a(y))i , which implies Π(x)i ≥ Π(y)i since Π(x)i = exp(ϕ(a(x)) i) P and Π(y)i = exp(ϕ(a(y))i ). On the other hand, from the definition of Π, we have Π(x)i = P Π(y)i = 1. Therefore, ϕ(a(x))i = ϕ(a(y))i for every i, and hence ψ(a(x))i − ψ(a(y))i = a(y)i − a(x)i ∈ [−M, −m]. Case 3: m ≤ M < 0. Similar to Case 2, we have ϕ(a(x))i ≤ ϕ(a(y))i for every i, and hence ψ(a(x))i − ψ(a(y))i = a(y)i − a(x)i ∈ [−M, −m]. Combining the three cases completes the proof. Lemma L.2 (Birkhoff contraction [Birkhoff, 1957]). Let A ∈ Rn×m ++ . Then dH (Ax, Ay) ≤ τ (A) dH (x, y), where τ (A) = tanh
∆(A) 4
< 1,
∆(A) = log max
Aik Ajl
i,j,k,l Ail Ajk
.
Lemma L.3 (Contraction of the composite map). Define Π K ⊤ (µ/(Kv)) . T (v) := K ⊤ (µ/(Kv)) v Then for all v, w ∈ Rk++ ,
dH (T (v), T (w)) ≤ τ (K ⊤ )τ (K) dH (v, w). v Proof. Let v, w ∈ Rk++ and define u(v) := µ/(Kv). Since log u(v) = log µ − log(Kv), we have
dH (u(v), u(w)) = dH (Kv, Kw). By the Birkhoff contraction lemma, dH (Kv, Kw) ≤ τ (K)dH (v, w),
dH (K ⊤ u(v), K ⊤ u(w)) ≤ τ (K ⊤ )dH (u(v), u(w)).
By Lemma L.1, dH
Π(K ⊤ u(v)) Π(K ⊤ u(w)) , K ⊤ u(v) K ⊤ u(w)
≤ dH (K ⊤ u(v), K ⊤ u(w)).
Combining the three inequalities yields the claim. Lemma L.4 (Hilbert-Lipschitz continuity of Π). Under the assumptions of Lemma L.1, for all v x, y ∈ Rk++ , dH (Π(x), Π(y)) ≤ 2dH (x, y). Proof. Let F (x) := Π(x)/x as in Lemma L.1, so Π(x) = x ⊙ F (x) componentwise. For a, b, c, d ∈ v Rk++ , a i bi ai bi ≤ max max , max i i i ci di ci di a i bi ai bi min ≥ min min . i ci di i ci i di Hence dH (a ⊙ b, c ⊙ d) ≤ dH (a, c) + dH (b, d). 28
Applying this with (a, b, c, d) = (x, F (x), y, F (y)) and using Lemma L.1, dH (Π(x), Π(y)) = dH (x ⊙ F (x), y ⊙ F (y)) ≤ dH (x, y) + dH (F (x), F (y)) ≤ 2dH (x, y).
Proof of Theorem 4.3. By Lemma L.3, T is a contraction in Hilbert projective metric with constant c := τ (K ⊤ )τ (K) < 1. Hence T has a unique fixed point in projective space (equivalently, unique up to positive scaling), v and for any initial v (0) ∈ Rk++ , dH (v (t) , v ∗ ) ≤ ct dH (v (0) , v ∗ ). This proves linear convergence of {v (t) }. It remains to show the KL bound for the output distributions. Let µ s(v) := K ⊤ . Kv By definition of T , for any fixed point v ∗ we have v ∗ = T (v ∗ ) =
Π(s∗ ) , s∗
s∗ := s(v ∗ ),
and with q ∗ := Π(s∗ ) this gives v ∗ = q ∗ /s∗ componentwise. Using dH (µ/(Kv), µ/(Kw)) = dH (Kv, Kw) and Birkhoff contraction twice, µ µ dH (s(t) , s∗ ) = dH K ⊤ , K⊤ (t) Kv ∗ Kv (t) ∗ t+1 ≤ c dH (v , v ) ≤ c dH (v (0) , v ∗ ). Since q (t) = Π(s(t) ), Lemma L.4 yields dH (q (t) , q ∗ ) = dH (Π(s(t) ), Π(s∗ )) ≤ 2dH (s(t) , s∗ ) ≤ 2ct+1 dH (v (0) , v ∗ ). Finally, because mj > 0 and q (t) , q ∗ ∈ Qm,ε , all components of q (t) , q ∗ are strictly positive. Let (t) ri := qi /qi∗ . Since both vectors are probability distributions, min ri ≤ 1 ≤ max ri . i
i
Therefore, KL(q
(t)
X (t) maxi ri = dH (q (t) , q ∗ ), ∥q ) = qi log ri ≤ log max ri ≤ log i min r i i i ∗
and the KL bound follows by combining with the previous inequality.
M
Additional Experiments
M.1
Experiments on Geometric Data
Here, we present additional experimental results on Gowalla location dataset [Cho et al., 2011]. In the Gowalla dataset, each user has a set of check-in locations, which we regard as a private distribution. We consider a scenario where each user privately releases a check-in location and a central server aggregates the released locations to estimate the overall distribution of check-ins. We use check-in data of 500 users in Austin, Texas and partitioned the area into a 20 × 20 grid and adopted the 29
Figure 4: Estimated distribution of check-ins for ε = 4. The true Figure 3: Comparison of distribution is shown in the leftmost panel, and the distributions WPM and KPM on Gowalla estimated by WPM and KPM are shown in the middle and right dataset. panels, respectively. Euclidean distance as a metric. Then, we evaluate the performance by the Wasserstein distance (p = 1) between the true distribution of check-ins and the estimated distribution from the released locations. For WPM, we used, as the base measure, the one that minimizes the worst-case cost among uniform measures. We report the mean and standard deviation over 5 runs with different choice of check-ins. Fig. 3 illustrates the comparison of the proposed WPM with KPM. We see that WPM consistently outperforms KPM across different privacy budgets. In addition, we also provide examples of the estimated distribution of check-ins for ε = 4 in Fig. 4. We observe that WPM captures the true distribution more accurately than KPM, which demonstrates the advantage of WPM in preserving the geometric structure of the data.
N
Experimental Details
Datasets Movielens-100K [Harper and Konstan, 2015] contains 100,000 ratings from 1000 users on 1700 movies, which is distributed under a usage license that permits research use. We randomly divided the movies into 80% training and 20% test sets: the training movies are used to construct each user’s input distribution (proportional to ratings), while the test movies form the candidate set for recommendation. Users are split into 80% validation and 20% test sets; the validation users are used for hyperparameter tuning, and the test users for final evaluation. Gowalla dataset [Cho et al., 2011] contains 6,442,892 check-ins from 196,591 users at 1,280,969 locations, which is distributed for academic research use. We used the check-in data in Austin, Texas. Computational environment Our experiments were conducted on Intel(R) Xeon(R) Silver 4316 CPU @ 2.30GHz and 8 NVIDIA A100-SXM4-80GB GPUs. Detailed setup for MovieLens experiments We computed the embeddings of movies using the top 64 principal components of the user-movie rating matrix on validation users. We use top 250 movies (based on the validation set of users) as the candidate set for recommendation. The privacy budget ε is swept over {0.2, 0.5, 1.0, 2.0, 3.0, 4.0}. The scale of the base measure for WPM and KPM is swept over the geometric grid in the feasible range [e−ε/2 /k, eε/2 /k]. We set λ = 0.01 and the number of iterations to 40 for WPM, which is sufficient for convergence in our experiments. Detailed setup for Gowalla experiments We use check-in data in Austin, Texas and drop users with fewer than 10 check-ins. For WPM, we use the optimal base measure among uniform measures and set λ = 0.01 and the number of iterations to 40, which is sufficient for convergence in our experiments.
O
Broader Impacts
This work aims to develop sampling mechanisms with formal privacy guarantees, and thus has the potential to enhance privacy in various applications. However, as with any privacy-preserving technology, there is a risk of misuse or unintended consequences. For example, if the privacy 30
parameters are not chosen carefully, it could lead to either insufficient privacy protection or excessive noise that degrades utility.
31