Locally Optimal Private Sampling: Beyond the Global Minimax

We thank the reviewer for their constructive feedback. We are glad they found the paper well-organized and appreciated our contribution to private sampling via a local minimax-optimal framework that extends beyond worst-case analysis. Below, we address each of the reviewer’s comments in detail.

1. The Condition $\\boldsymbol{\frac{c_2 - c_1}{1 - c_1} \in \mathbb{N}}$:

We define $\tilde{\mathcal{P}}$ as the set of distributions of interest over a sample space $\mathcal{X}$. To enable a unified theoretical treatment of both discrete and continuous settings, we restrict $\tilde{\mathcal{P}}$ to be $\tilde{\mathcal{P}}_{c_1,c_2,\mu}:=\\{P\in\mathcal{P}(\mathcal{X}): P \ll \mu,\\; c_1\leq\frac{dP}{d\mu}\leq c_2\\; \mu\text{-a.e.}\\}$ for some measure $\mu$. As shown in Park et al. [17] and mentioned in the manuscript, such restriction is necessary to avoid the trivial result. For this general setting, we assume the technical condition $\frac{c_2-c_1}{1-c_1} \in \mathbb{N}$ to facilitate the construction of our proofs via a decomposability argument.

In the discrete case, where the domain is $[k]$, we define $\mathcal{P}([k])=\tilde{\mathcal{P}}_{0,k,\mu_k}$, where $\mu_k$ is the uniform distribution over $[k]$. In this setting, we have $\mathcal{X}=[k]$, $c_1=0$, $c_2=k$, and $\mu=\mu_k$, which yields $\frac{c_2-c_1}{1-c_1}=k \in\mathbb{N}$. Thus, the above assumption trivially holds.

In the continuous case, following the approach of Park et al. [17], we assume $\tilde{\mathcal{P}} \subset \tilde{\mathcal{P}}\_{c_1,c_2,\mu}$ with base measure $\mu$ being Lebesgue. Here, $\tilde{\mathcal{P}}_{c_1, c_2, \mu}$ acts as a tractable superset of the true distribution class. If $\frac{c_2-c_1}{1-c_1}\notin\mathbb{N}$ for the original parameters, one can slightly increase $c_2$ or decrease $c_1$ to satisfy the condition. This adjustment results in a slightly enlarged superset, which is still consistent with the modeling assumption that $\tilde{\mathcal{P}}$ lies within a known, manageable class of distributions.

For example, setting $c_1=0$ always ensures the condition holds for any integer $c_2$. Moreover, in some natural families of distributions—such as mixtures of Gaussians with shared variance (see Example 3.1)—the assumption is satisfied without needing any relaxation.

In summary, the condition $\frac{c_2-c_1}{1-c_1}\in\mathbb{N}$ is not restrictive: it is either automatically satisfied in discrete settings, or can be enforced without loss of generality in continuous ones, purely to streamline the technical development of the theory.

2. Intuition Behind the Outperformance of the Non-Linear Sampler:

We thank the reviewer for the insightful suggestion. The non-linear sampler outperforms the linear one because it is not only worst-case optimal but also instance-optimal: for each input distribution $P \in \tilde{\mathcal{P}}_{c_1,c_2,\mu}$ and any $f$-divergence $D_f$, it minimizes $D_f(P\\|\mathbf{Q})$ over all admissible $\mathbf{Q}$ satisfying $\varepsilon$-LDP (see Proposition D.2, line 1273). This instance-optimality arises because the non-linear sampler is constructed by solving the optimization problem for each individual $P$, producing a privatized distribution that minimizes the divergence from $P$ under the given LDP constraints. In contrast, the linear sampler applies the same transformation across all inputs. To provide clearer context for the reader, we will incorporate this intuition before presenting the comparison between the linear and non-linear samplers in the main manuscript.

3. Practical Utility of the Proposed Samplers and Discussion of Potential Applications:

We thank the reviewer for the suggestion. We agree that including concrete downstream applications helps underscore the potential impact of our framework.

Our framework is particularly relevant in settings where private samples must be generated in the presence of public data. For example, recent work on differentially private next-token prediction (e.g., Flemings et al. [28]) forms sampling distributions by linearly combining private and public distributions (e.g., from a public LLM), closely resembling our linear samplers. Our approach offers a principled, theoretically grounded alternative with formal optimality guarantees.

As discussed in the introduction, our work is also motivated by the growing relevance of personalized collaborative learning [18–20], where individual models are trained collaboratively and public data may reflect privatized information shared by other users. Since real-world distributions are rarely worst-case, our local minimax formulation provides a more realistic and effective framework for capturing achievable performance in practice. This makes our approach applicable to tasks such as private model fine-tuning, federated estimation, and privacy-preserving generative modeling.

We acknowledge that our experiments are limited to synthetic data due to the computational overhead of the clipping step in high-dimensional settings, as noted in the limitations section (lines 328–333). A promising direction for future work is to develop a scalable approximation for the optimal sampler that does not require clipping. Such approximation would enable us to implement private sampling on high-dimensional real datasets.

4. Connection to Existing Literature in LDP Estimation:

We thank the reviewer for pointing this out. Several prior works have explored local minimax formulations under DP, primarily in estimation settings. For example, Chen et al. [45] study discrete distribution estimation under local DP using neighborhoods around a reference distribution. McMillan et al. [47] and Duchi et al. [43] develop locally minimax estimators for 1D parameters in the central and local models, respectively. Asi et al. [48] extend this to function estimation, and Feldman et al. [44] study nonparametric density estimation with both local and instance-optimal guarantees.

Thus, local minimax formulations using neighborhood structures have been well-studied in DP estimation and distribution learning, with definitions tailored to specific inference goals. As outlined in the related work section, while our focus differs in that we study private sampling rather than estimation, our formulation adopts a similar neighborhood-based perspective, adapted to the sampling setting. We will revise the related work section to clarify this connection and include relevant context.

5. Computational Complexity of Clipping in High Dimensions:

The optimal sampler in Theorem 5.1 is defined by the mapping $\mathbf{Q}^\star_\varepsilon(P)$, whose density is

where $r_P>0$ is chosen so that $\int_{\mathbb{R}^n}q(x)dx=1$. Therefore, the computational burden of $\mathbf{Q}^\star_\varepsilon(P)$ is dominated by computing $r_P$. To that end, consider the family $q_r(x):=\text{clip}\bigl(\tfrac{1}{r}p(x);\tfrac{\gamma+1}{\gamma+e^{\varepsilon}}p_0(x),\tfrac{(\gamma+1)e^{\varepsilon}}{\gamma+e^{\varepsilon}}p_0(x)\bigr)$ and search for $r>0$ satisfying $\int_{\mathbb{R}^n}q_r(x)dx=1$. Let $r_1:=\tfrac{\gamma+e^{\varepsilon}}{\gamma(\gamma+1)}$ and $r_2:=\tfrac{\gamma(\gamma+e^{\varepsilon})}{e^{\varepsilon}(\gamma+1)}$. It can be shown that $\int_{\mathbb{R}^n}q_{r_1}(x)dx\ge1$ and $\int_{\mathbb{R}^n}q_{r_2}(x)dx\le1$, while the map $r\mapsto\int q_r(x)dx$ is continuous. Therefore, the desired normalizing constant $r_P$ can be computed via the bisection method with initial endpoints $(r_1,r_2)$.

While the bisection procedure is efficient, the repeated numerical integration required at each step can become computationally prohibitive in high-dimensional settings. This underscores the need for alternative approaches that avoid explicit integration. One promising direction is to design approximate samplers by minimizing $D_f(\mathbf{Q}\\|\mathbf{Q}^\star_\varepsilon(P))$, where $\mathbf{Q}^\star_\varepsilon$ is the optimal $\varepsilon$-LDP sampler.

We consider these directions to be promising avenues for future work, as they address an important and nontrivial challenge in scaling optimal clipping-based mechanisms to high-dimensional domains.

6. Connection to Instance-Optimal Differential Privacy:

The concept of instance-optimality, formulated by Feldman et al. [44], refers to algorithms that perform nearly as well on a given instance as the best differentially private algorithm designed specifically for that instance and its neighborhood. To avoid trivial solutions, performance is benchmarked not just against the instance itself but over a surrounding neighborhood, ensuring meaningful generalization without requiring the algorithm to be tailored to each instance.

Instane-optimality captures the idea that algorithms should adapt to the "hardness" of the input rather than the worst case—particularly useful when there is a large gap between typical and worst-case performance, as in nonparametric estimation.

In contrast, our local minimax approach evaluates worst-case risk over a neighborhood around a reference distribution, offering robust guarantees across plausible inputs rather than individual adaptation. This is especially relevant when prior knowledge or public data (e.g., a shared $P_0$) is available to guide the analysis.

That said, although the two formulations address different goals, they share structural similarities. As noted in the related work section, the neighborhood $N(P_0)$ in our framework aligns with that in Feldman et al.’s instance-optimality. We will add this comparison to the appendix to clarify their conceptual and technical connections.

Thank you so much for clarifying these issues. While I appreciate the efforts and may consider increasing my score, I remain concerned about the practical relevance of the privacy guarantee.

Thank you for the clarification.

We thank the reviewer for their thoughtful review and insightful question. We are glad that they found the paper well written and appreciated our contributions, including the generalization of the global minimax framework and the introduction of instance-optimal samplers. Below, we provide a detailed discussion of the distribution class used in our work.

The class of distributions for which the sampler is optimal is restricted (the density verifies: $c_1 h(x) \leq p(x) \leq c_2 h(x)$). Although some restrictions are necessary to exclude trivial samplers, it is unclear to me why this particular condition arises or whether other choices could have been made.

Could you discuss the choice of distribution class?

We thank the reviewer for raising this important point. The condition $c_1 h(x) \leq p(x) \leq c_2 h(x)$ used to define the distribution class $\tilde{\mathcal{P}}_{c_1, c_2, h}$ may seem restrictive at first glance, but it arises naturally from the need to balance generality with analytical tractability in LDP sampling.

As discussed in Park et al. [17], one may consider broader classes such as: (i) $\mathcal{P}(\mathbb{R}^n)$, the set of all probability measures on $\mathbb{R}^n,$ or (ii) $\mathcal{C}(\mathbb{R}^n)$, the set of all continuous distributions on $\mathbb{R}^n$. However, for such general choices, it has been shown that the minimax rate is vacuous, that is, any $\varepsilon$-LDP sampling mechanism attains the largest possible value of $f$-divergence in the worst case (see Proposition 3.2 in Park et al.). As a result, all mechanisms perform equally poorly under these settings, making meaningful comparison or claims of optimality infeasible.

To obtain nontrivial and practically useful bounds, one must restrict attention to more regular classes of distributions. The class $\tilde{\mathcal{P}}_{c_1, c_2, h}$, which consists of distributions whose densities (with respect to a general base measure such as Lebesgue) are bounded between two scalings of a reference function $h$, is one such choice. This class offers two key properties: (i) it includes a wide range of realistic distribution families (so it is rich enough), as we illustrate below, and (ii) it ensures that the worst-case $f$-divergence is non-vacuous.

This modeling choice is not purely theoretical. As shown in Example 3.1 of our manuscript, mixtures of Gaussian or Laplace distributions with bounded location parameters fall within this class. For instance, Laplace mixtures supported within the $\ell_1$ unit ball satisfy $\\tilde{\mathcal{P}}_{\\mathcal{L}} \\subseteq \\tilde{\\mathcal{P}}\_{e^{-1/b}, e^{1/b}, h\_{\\mathcal{L}}}$ and Gaussian mixtures with bounded means lie in $\\tilde{\\mathcal{P}}\_N \subseteq \\tilde{\\mathcal{P}}\_{0, 1, h\_{\\mathcal{N}}},$ for suitable choices of the reference functions $h\_{\\mathcal{L}}$ and $h\_{\\mathcal{N}}$ as defined in Example 3.1. These examples demonstrate that the proposed model captures a broad and practically relevant class of distributions.

In summary, while other distribution classes are conceivable, the choice of $\tilde{\mathcal{P}}_{c_1, c_2, h}$ provides a balanced trade-off between expressiveness and analytical feasibility, and encompasses many distribution families relevant to real-world applications.

We thank the reviewer for their thoughtful feedback and for pointing out areas needing clarification, such as the definitions of $g^*$ and $\tilde{\mathcal{P}}$. We are also glad they found the extension from global to local minimax meaningful and acknowledged our contribution to private distribution sampling. Below, we address each of the reviewer’s questions and comments individually.

1. Presentation and Organization of Section 3:

We appreciate the reviewer’s suggestion regarding the presentation and structure of Section 3. While the main focus of the paper is on local sampling, we included Section 3 to provide the necessary theoretical groundwork for the subsequent sections, including the extension to functional LDP, the formulation of key assumptions, and the definition of global samplers used in the analysis. Specifically, we study global minimax samplers under the general framework of functional LDP and derive compact characterizations of the optimal samplers. As noted at the beginning of Section 3 (lines 170–171), these global samplers serve as proxies for constructing local minimax-optimal samplers in the following sections.

We aimed to preserve a natural progression in the paper’s structure, moving from functional LDP global minimax in Section 3, to functional LDP local minimax in Section 4, and then to pure LDP local minimax in Section 5. This order reflects the way the results were derived, with each step building on the previous one. That said, we agree that the purpose of Section 3 could be emphasized more clearly. In the revision, we will improve the framing of this section and better articulate its connection to the rest of the paper. We will also consider reorganizing key parts into a dedicated subsection to enhance clarity and flow.

2. Computational Efficiency of the Proposed Sampler and Clarification of the Projection Step in Theorem 4.1:

For general $f$-divergences and sample spaces $\mathcal{X}$, this projection step may require solving a constrained optimization problem, and its complexity is on the same order as that of the clipping-based projection described in Theorem 5.1. However, as shown by Husain et al. [7], in certain special cases such as when the $f$-divergence is the KL divergence and the sample space is finite, the projection problem admits a closed-form solution using the Karush-Kuhn-Tucker (KKT) conditions. In these settings, the projection step can be carried out efficiently without incurring additional computational overhead.

As the reviewer points out, the projection can be computationally demanding, particularly in high-dimensional settings, even though a closed-form expression for it is available. This is because the expression involves a clipping step that depends on a normalization constant (similar to $r_P$ in the statement of Theorem 5.1). Computing this constant requires performing numerical integration over a high-dimensional space, which can be costly and often becomes the main computational bottleneck.

In the following, we include a discussion on the computational complexity of the samplers proposed in our work. Overall, the samplers derived in this work fall into two categories: (1) linear samplers, such as the core sampler in Theorem 4.1, and (2) nonlinear samplers involving clipping functions, as in Theorem 5.1.

The linear sampler $\mathbf{Q}\_{c_1,c_2,h,g}^\star$, defined in Theorem 3.4 (Equation (4)), has density $q_g^{\star}(P)(x) := \lambda_{c_1, c_2, g}^{\star} p(x) + \left(1 - \lambda_{c_1, c_2, g}^{\star} \right) h(x)$, a convex combination of the input $p$ and reference $h$. Its complexity lies in computing the mixing coefficient $\lambda_{c_1,c_2,g}^{\star} = \inf_{\beta\geq0}\frac{e^{\beta}+\frac{c_2-c_1}{1-c_1}\left(1+g^*(-e^{\beta})\right)-1}{(1-c_1)e^{\beta}+c_2-1}$.

In many cases—e.g., pure, approximate, and Gaussian LDP—the conjugate $g^*$ is known in closed form, allowing this optimization to be solved analytically or efficiently numerically. Thus, the linear sampler is computationally tractable.

In contrast, the nonlinear sampler in Theorem 5.1 is defined by the mapping $\mathbf{Q}^\star_\varepsilon(P)$, whose density is

where the constant $r_P>0$ is selected to ensure normalization condition $\int_{\mathbb{R}^n} q(x) dx = 1$. The main computational burden lies in determining this normalizing constant $r_P$. To compute it, we define the parametric family $q_r(x) := \operatorname{clip}\left( \tfrac{1}{r}p(x); \tfrac{\gamma+1}{\gamma+e^{\varepsilon}}p_0(x), \tfrac{(\gamma+1)e^{\varepsilon}}{\gamma+e^{\varepsilon}}p_0(x)\right)$, and search for a value $r > 0$ such that $\int_{\mathbb{R}^n} q_r(x) dx = 1$. It can be shown that the endpoints $r_1 := \tfrac{\gamma+e^{\varepsilon}}{\gamma(\gamma+1)}$ and $r_2 :=\tfrac{\gamma(\gamma+e^{\varepsilon})}{e^{\varepsilon}(\gamma+1)}$ satisfy $\int q_{r_1}\ge 1$ and $\int q_{r_2} \le 1$, and that $r \mapsto \int q_r(x)\,dx$ is continuous. Therefore, $r_P$ can be found using the bisection method on the interval $(r_1,r_2)$.

While the bisection method requires only a logarithmic number of steps, the main computational cost lies in evaluating the integral $I(r) = \int_{\mathbb{R}^n} q_r(x) dx$ at each step. This is tractable in low-dimensional settings but quickly becomes prohibitive as dimensionality increases.

Regarding the reviewer’s question about the expression “$\hat{P}$ minimizes $\inf_{P' \in N(P_0)} D_f(P \\| P')$,” we intended to convey that we find the distribution $\hat{P}$ that minimizes $D_f(P \\| P')$ over all $P' \in N(P_0)$. The reviewer is correct, and we will revise the text accordingly for clarity.

3. Relationship Between Global and Local Sampling and the Associated Rate Gap:

Under an $\varepsilon$-LDP constraint, consider the global model class

where $r_1 = \frac{e^{\varepsilon} + \gamma}{\gamma(\gamma + 1)}$ and $r_2 = \frac{\gamma(e^{\varepsilon} + \gamma)}{e^{\varepsilon}(\gamma + 1)}.$

While the difference between the local and global risks can be computed in closed form via direct subtraction, the resulting expression is not particularly insightful, as it depends on several problem-specific parameters, including $(c_1, c_2, \gamma)$, the sample space, and the privacy level $\varepsilon$. Although the subtraction can be evaluated for any given instance, it does not, in general, yield a simple or interpretable expression.

To better illustrate the improvement, we provide numerical comparisons in both discrete and continuous settings. In the discrete case (Figures 3, 5, 6), we compare the global and local minimax-optimal samplers for $\mathcal{X} = [k]$, where the global class is $\tilde{\mathcal{P}}\_{\text{global}} = \mathcal{P}([k]) = \tilde{\mathcal{P}}\_{0,k,\mu_k},$ with $\mu_k$ being the uniform distribution. The local neighbourhood is defined as $\tilde{\mathcal{P}}\_{\text{local}} = \mathcal{N}\_\gamma(\mu_k) = \tilde{\mathcal{P}}\_{\frac{1}{\gamma}, \gamma, \mu_k},$ with $\gamma = \frac{k}{2} - 1$, ensuring that $\mathcal{N}_\gamma(\mu_k) \subseteq \mathcal{P}([k])$.

We instantiate both the global and local neighbourhoods using the parameters defined in Section 6.1, and compute the exact minimax risks in each case. Figures 3, 5, and 6 present these theoretical worst-case risks, clearly demonstrating the benefit of tailoring the mechanism to the local neighbourhood. These plots consistently illustrate a significantly improved minimax rate in the local setting compared to the global one, across various $f$-divergences and parameter ranges.

4. Definition of $\boldsymbol{g^*}$ in Equation (4):

We define $g^{\*}$ to be the convex conjugate of the function $g$, as introduced in the notation section. Formally, it is given by $g^{*}(y) = \sup_{x \in \text{dom}(g)} \Bigl( xy - g(x) \Bigr)$. To improve clarity and completeness, we will include this formal definition in the notation section.

5. Definition of $\boldsymbol{\tilde{\mathcal{P}}}$ in Section 4:

As for the definition of $\tilde{\mathcal{P}}$ in Section 4, we clarify that $\tilde{\mathcal{P}}$ refers to a general set of possible distributions relevant to the problem under consideration. In particular, in Section 4, $\tilde{\mathcal{P}}$ may represent any general class of continuous distributions, subject only to the constraint $N_\gamma(P_0) \subseteq \tilde{\mathcal{P}},$ which is assumed in the statement of Theorem 4.1. We will include a note at the beginning of the section to make this generality of $\tilde{\mathcal{P}}$ explicit.

We thank the reviewer for their thoughtful and detailed review. We are glad they appreciated our contributions, including the extension to functional LDP and the characterization of optimal local minimax samplers. Below, we address each of the reviewer’s questions and concerns individually.

The focus of the paper on the local model makes the sampling problem less about algorithms and more about finding the right parameters.

While global minimax private sampling under pure LDP has been studied previously (e.g., Park et al. [17]), our work extends this by developing both global and local minimax results under the more general functional LDP framework, which encompasses pure, approximate, and Gaussian LDP. This required introducing new assumptions, proof techniques, and a carefully defined local neighborhood that fundamentally differs from the global setting. These developments allow us to bridge global and local minimax frameworks and derive new linear and nonlinear samplers that are optimal in the local minimax sense under functional LDP. Importantly, we propose new algorithms tailored to this setting—not merely new parameters for existing ones. Our framework also lays a foundation for future work, including broader neighborhood definitions and scalable implementations, as outlined in the limitations.

The paper doesn't give much intuition about the derived samplers.

We thank the reviewer for their constructive feedback. While we made an effort to convey the intuition behind the derived samplers—through both textual explanations and visual illustrations—we agree that additional elaboration could help improve clarity. Below, we provide a complementary interpretation of the two main classes of samplers discussed in the paper: linear and nonlinear samplers.

For the linear samplers, such as the global sampler $\mathbf{Q}^\star_{c_1,c_2,h,g}$ defined by $q^\star_g(P)(x):=\lambda^\star_{c_1,c_2,g}p(x)+\left(1-\lambda^\star_{c_1,c_2,g}\right)h(x),$ as presented in lines 204–206 of the manuscript, a natural interpretation is that the mechanism samples from the input distribution $P$ with probability $\lambda^\star_{c_1,c_2,g}$, and from the reference distribution $h$ otherwise. This probabilistic mixture reflects the privacy–utility trade-off: a larger $\lambda^\star_{c_1,c_2,g}$ yields outputs closer to $P$, while sampling from $h$ introduces the randomness required for privacy. This same intuition applies to all linear samplers derived in the paper.

For the nonlinear sampler presented in Theorem 5.1, the optimal sampling mechanism is given by
$\hspace{6cm}q(x)=\text{clip}\left(\tfrac{1}{r_P}p(x);\tfrac{\gamma+1}{\gamma+e^{\varepsilon}}p_0(x), \tfrac{(\gamma+1)e^{\varepsilon}}{\gamma+e^{\varepsilon}}p_0(x)\right).$

As described in lines 287–288 of the manuscript, this construction defines a nonlinear sampler by projecting the input distribution $P$ onto a family of distributions referred to as the mollifier $M_\varepsilon$, as defined in the paper (see Figure 2). The clipping operation ensures that the output distribution $\mathbf{Q}^\star_\varepsilon(P)$ satisfies the LDP constraint while maintaining a close match to the original distribution $P$ in terms of $f$-divergence. In other words, $\mathbf{Q}^\star_\varepsilon(P)$ belongs to the mollifier and is instance-optimal for the given input distribution. This approach yields a nonlinear transformation of $P$ that carefully balances privacy constraints and utility preservation.

Figure 2 further illustrates the sampling process under the local minimax framework. It involves two steps: (1) if $P$ lies outside the local neighborhood of the public distribution $P_0$, it is projected inward; and (2) optimal sampling is performed using either a linear or nonlinear sampler, depending on the context.

We will incorporate these clarifications and interpretations in the final version to better communicate the intuition behind the derived samplers.

Can you please provide a comparison between how you use public data and previous work [22] used them?

Zamanlooy et al. [22] use public data as an explicit optimization constraint in their global minimax formulation, which is defined over discrete domains and restricted to linear samplers. Note that in this discrete domain, the sampling distribution can be expressed as $PK$, where $P$ is the data distribution (viewed as a vector) and $K$ is an LDP mechanism (viewed as a matrix that satisfies the LDP constraint). Zamanlooy et al. [22] considered $\varepsilon$-LDP mechanisms $K$ that preserve the public prior $P_0$, enforcing $P_0K = P_0$. This guarantees that the public data remains invariant under the mechanism while privacy is provided for private data. They then formulated their minimax optimization as: $\hspace{8cm}\inf_{\\substack{K \in\mathcal{Q_\varepsilon}\\\\P_0K=P_0}}\sup_{P \in\mathcal{P}(\mathcal{X})}D_f(P\\|PK).$

In contrast, we use public data in a conceptually different way. Instead of treating it as a hard constraint, we use the public distribution $P_0$ to define a local neighborhood $N_\gamma(P_0)$. The worst-case input distributions are then considered within this neighborhood. This leads to a local minimax formulation of the private sampling problem over a general sample space, without restricting to linear mechanisms. Specifically, we search over mechanisms $\mathbf{Q}\in\mathcal{Q_\varepsilon}$, where $\mathcal{Q_\varepsilon}$ again denotes the class of $\varepsilon$-LDP samplers:

Note that the formulation of [22] is a special form of our setting: $\mathbf{Q}(P)=PK$ where $K$ is an $\varepsilon$-LDP sampler. Our goal is to move beyond the overly pessimistic global worst-case risk and instead focus on distributions that are close to $P_0$, which is assumed to be publicly available to all users. This enables a more refined and robust characterization of optimal samplers from a local minimax perspective. It is important to note that in the formulation of Zamanlooy et al., having public data does NOT necessarily lead to a more accurate sampling (smaller minimax risk), whereas this is the case in our formulation.

That said, as shown in Theorems 4.1 and 5.1, the local minimax-optimal samplers $Q^\star_{g, N_\gamma(P_0)}$ and $Q^\star_{\varepsilon,N_\gamma(P_0)}$ both satisfy $Q^\star_{g, N_\gamma(P_0)}(P_0)=P_0$ and $Q^\star_{\varepsilon, N_\gamma(P_0)}(P_0)=P_0$, even though the preservation of $P_0$ is not explicitly enforced as a constraint. This demonstrates that our formulation inherently respects the public data when the input distribution coincides with it, without requiring this condition to be imposed as part of the optimization, meaning that our formulation is a non-trivial generalization of Zamanlooy et al. [22] to general data distribution (not necessarily to discrete distribuition, as in [22]).

In the limitation you mentions the computational complexity of your samplers, can you give more details on it?

We thank the reviewer for raising this important point regarding computational complexity. Overall, the samplers derived in this work fall into two categories: (1) inear samplers, such as the core sampler in Theorem 4.1, and (2) nonlinear samplers involving clipping functions, as in Theorem 5.1.

The linear sampler \\mathbf{Q}^\\star_{c_1,c_2,h,g} defined in Theorem 3.4 (Equation (4)), has density q^\\star_g(P)(x):=\\lambda_{c_1,c_2,g}^{\\star}\\,p(x)+\\left(1-\lambda_{c_1,c_2,g}^{\star} \right) h(x), a convex combination of the input $p$ and reference $h$. Its complexity lies in computing the mixing coefficient $\lambda_{c_1,c_2,g}^{\star}$, which is given by $\lambda_{c_1,c_2,g}^{\star} = \inf_{\beta\geq0}\frac{e^{\beta}+\frac{c_2-c_1}{1-c_1}\left(1+g^*(-e^{\beta})\right)-1}{(1-c_1)e^{\beta}+c_2-1}$.

In many cases—e.g., pure, approximate, and Gaussian LDP—the conjugate $g^*$ is known in closed form, allowing this optimization to be solved analytically or efficiently numerically. Thus, the linear sampler is computationally tractable.

In contrast, the nonlinear sampler in Theorem 5.1 is defined by the mapping $\mathbf{Q}^\star_\varepsilon(P)$, whose density is $\hspace{6cm}q(x)=\text{clip}\left(\tfrac{1}{r_P}p(x);\tfrac{\gamma+1}{\gamma+e^{\varepsilon}}p_0(x), \tfrac{(\gamma+1)e^{\varepsilon}}{\gamma+e^{\varepsilon}}p_0(x)\right),$ where the constant $r_P>0$ is selected so that $\int_{\mathbb{R}^n}q(x)dx=1$. The main computational burden lies in determining this normalizing constant $r_P$. To compute it, we define the parametric family $q_r(x):=\text{clip}\left( \tfrac{1}{r}p(x); \tfrac{\gamma+1}{\gamma+e^{\varepsilon}}p_0(x),\tfrac{(\gamma+1)e^{\varepsilon}}{\gamma+e^{\varepsilon}}p_0(x) \right)$, and search for a value $r>0$ such that $\int_{\mathbb{R}^n}q_r(x)dx=1$. It can be shown that the endpoints $r_1:=\tfrac{\gamma+e^{\varepsilon}}{\gamma(\gamma+1)}$ and $r_2:=\tfrac{\gamma(\gamma+e^{\varepsilon})}{e^{\varepsilon}(\gamma+1)}$ satisfy $\int q_{r_1}\ge1$ and $\int q_{r_2} \le 1$, and that $r \mapsto \int q_r(x)dx$ is continuous. Therefore, $r_P$ can be found using the bisection method on the interval $(r_1,r_2)$.

While the bisection method requires only a logarithmic number of steps, the main computational cost lies in evaluating the integral $I(r)=\int_{\mathbb{R}^n}q_r(x)dx$ at each step. This is tractable in low-dimensional settings but quickly becomes prohibitive as dimensionality increases.

A promising direction to mitigate this is to bypass explicit integration by designing approximate samplers $\mathbf{Q}$ that minimize the divergence $D_f(\mathbf{Q}\\|\mathbf{Q}^\star_\varepsilon(P))$, where $\mathbf{Q}^\star_\varepsilon$ is the optimal $\varepsilon$-LDP sampler. This offers a rich path for future work, especially in scaling optimal mechanisms to high-dimensional domains.