Differentially Private Nonparametric Modal Learning with Applications to Regression and Clustering

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Differentially Private Nonparametric Modal Learning with Applications to Regression and Clustering

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

Differentially Private Nonparametric Modal Learning with Applications to Regression and Clustering Arkajyoti Bhattacharjee and Arnab Auddy Department of Statistics The Ohio State University August 3, 2026

arXiv:2607.29675v1 [math.ST] 31 Jul 2026

Abstract Density modes provide a localized and interpretable summary of complex, multimodal distributions, but their estimation under rigorous differential privacy constraints remains largely unexplored. We study differentially private recovery of density modes for multivariate distributions under local smoothness, curvature, and separation conditions. We propose DP-GRAMS, a mean-shift inspired method that performs noisy ascent on a differentially private score estimator. Assuming the density belongs locally to a Hölder class with smoothness parameter β > 2, our score estimator uses bias-reducing higher-order kernels, and then enforces privacy in the gradient ascent steps via gradient clipping and calibrated Gaussian noise. A private initialization scheme combines a density-aware utility with a diversity-inducing suppression rule and, with k ≍ M log n draws over a public hDAP -grid and suppression radius ρinit ≍ (log n)−1/d , achieves high-probability coverage of the modal basins by successively suppressing selected local neighborhoods in competitive regions, while correlated noise across multiple starts enables joint release under a single (ε, δ)-differential privacy guarantee. We prove that all population modes are recovered and establish asymptotic error rates of the form with high probability 2(β−1) β−1 d+β O ( logn n ) d+2β + O ( polylog(n,δ) ) . We also provide minimax lower bounds for private 2 2 n ε mode estimation, and show that our estimators are nearly optimal, up to a logarithmic factor in the MSE. We present two natural extensions: DP-PMS, a private modal-regression method, and DP-GRAMS-C, a clustering pipeline. Extensive experiments on synthetic and real data demonstrate favorable privacy–utility trade-offs relative to common baselines.

Introduction

Estimating the modes of a probability density – its local maximizers – is a central problem in statistics. Unlike global summaries such as the mean or median, modes reveal heterogeneous subpopulations and localized concentrations of probability mass (Chacón, 2015; Chen et al., 2016a), making them indispensable in multimodal or complex settings. Applications span a wide range of domains, from clustering and classification (Avidan, 2007; Chen et al., 2016b; Li et al., 2007) to computer vision tasks such as object tracking and image segmentation (Comaniciu and Meer, 2002; Comaniciu et al., 2003), as well as nonlinear statistical modeling paradigms including manifold 1

learning and modal regression (Einbeck and Tutz, 2006; Chen et al., 2016a). Recovering modes is inherently nonregular: small perturbations of the density may create or destroy critical points, and statistical difficulty depends delicately on local smoothness, curvature, and separation (Tsybakov, 2008; Genovese et al., 2014). In this paper, we study the problem of mode estimation while maintaining privacy of individual data points. This is motivated from a practical standpoint, where computing distribution summaries from sensitive data that routinely arise in domains such as healthcare or finance presents significant privacy challenges. Recent research has shown that even such summaries can compromise privacy (Oberski and Kreuter, 2020; Dick et al., 2023). Differential privacy (Dwork et al., 2006, 2014) provides rigorous protection by ensuring that the output of an algorithm is nearly indistinguishable with or without any single individual. Considerable research has developed differentially private methods for means, regression, and clustering (see Alabi et al., 2020; Cai et al., 2021; Oberski and Kreuter, 2020; Dankar and El Emam, 2013), but private mode estimation has received little direct attention, with Pacchiano et al. (2021) as a notable exception. This motivates the current work, where we pose private mode estimation as a gradient-ascent problem and develop practically implementable techniques with statistical guarantees. We use the natural characterization of modes as critical points of the score function ∇x log p(x), and identify modes among these critical points through local Hessian conditions. In order to achieve this, we use a suitable estimate of the density, denoted pb(x), and the score function induced by it. To decrease the sensitivity of our score estimator from extreme observations, rather than working directly with ∇x log pb(x), we use a stabilized score estimator, obtained by clipping the KDE gradient and flooring the KDE denominator. More precisely, for parameters A > 0 and pfloor > 0, we define clipA (∇pb(x)) , sbA,pfloor (x) = max{pb(x), pfloor }

A clipA (v) := v min 1, ∥v∥2

where

(1)

To find critical points of the log-density, given an initialization, we use the estimated gradients at each iterate to run the gradient ascent steps: xt+1 = xt + η sbA,pfloor (xt ), where η is the stepsize. Our density and score estimation are based on kernel smoothing. Given P i data X1 , . . . , Xn ∈ Rd and kernel K, we use the kernel density estimator pb(x) = nh1 d ni=1 K x−X h and its gradient ∇pb(x) to compute our score estimator. We note here that scaled gradient ascent on log pb for Gaussian kernels coincides with mean shift (Fukunaga and Hostetler, 1975; Cheng, 1995; Comaniciu and Meer, 2002). Since we allow for multiple modes, we assume that the log density is concave locally, but not necessarily globally. To tackle estimation in this setting, we use a multi-tiered algorithm. First, we pick several initializations that are likely to be candidates for the mode. Second, we refine these initializations via the gradient ascent procedure mentioned above. Finally, we merge sufficiently close refinements to their local centers, thus allowing for a unified mode estimation process. To ensure differential privacy, we use several different mechanisms in the various stages of the algorithm: 2

1. Initialization: we use a density-aware private initialization scheme called DAP, which computes a local empirical-mass utility over a public grid and samples anchor points from the grid via the exponential mechanism (McSherry and Talwar, 2007). In addition to being density aware due to the exponential mechanism being tuned to the empirical-mass utility, the algorithm uses a local suppression step after each selected anchor, so that the same neighborhood is not repeatedly selected and new high density regions are visited. With k ≍ M log n draws, this scheme yields high-probability coverage of the neighborhood of each mode. 2. private stochastic gradient ascent: we use the stabilized score estimator sbA,pfloor in (1) to run a stochastic gradient ascent and at each gradient ascent step, add suitably calibrated Gaussian noise. Similar ideas have been used earlier in parametric problems, via DP-SGD (Bassily et al., 2014). 3. correlated noise: Since iterate sequences from multiple initializations can be close to each other, our privacy preserving noise is correlated across initializations. At iteration t, the correlation structure is induced by an exponential kernel on the current iterate locations xt,1 , . . . , xt,k . This ensures that multiple initializations do not lead to undue privacy loss. We refer to the entire mechanism above as DP-GRAMS (Differentially Private Gradient Ascent for Mode Seeking). We ensure that it satisfies (ε, δ)-differential privacy (see, e.g., Dwork, 2006). The privacy budget ε is appropriately apportioned into the initialization and gradient ascent steps mentioned above. c contains, with high probability, On the utility side, we show that the final merged estimator M bj satisfying for each true population mode µj , a released point µ 2

b j − µj ∥ ≲ E∥µ

log n n

2(β−1) d+2β

polylog(n, δ) n 2 ε2

β−1 d+β

Our results depend crucially on appropriate curvature assumptions and β-Hölder smoothness of the density function. Most importantly, we assume that modes are sufficiently separated and that the density is locally log-concave around each mode. Next, we use higher order kernels that can leverage advantages due to β-Hölder smoothness for β > 2. We complement the above results with the minimax lower bound inf

sup

b x∈T (n,ε,δ) p∈Pβ (L), µ∈modes(p)

b−µ Ep x

2(β−1)

≳ n

− d+2β

+ (nε)

−

2(β−1) d+β

where δ = o(n−1 ) and T (n, ε, δ) is the set of all possible estimators based on a sample of size n and satisfying (ε, δ) differential privacy. Thus our estimators are nearly minimax optimal, since the MSE upper bounds match the lower bounds up to logarithmic terms. As expected, the MSE separates into a non-private rate given by the first term (matching Theorem 3 of Tsybakov, 1990) while the second captures the degradation required for privacy. The recovered private modes also motivate downstream procedures for modal regression and clustering. We develop: (i) DP-PMS, which adapts the private ascent mechanism to conditional 3

mode estimation by updating in the response direction, and (ii) DP-GRAMS-C, which releases private modal centers for clustering and uses deterministic assignments as post-processing or evaluation. These procedures use the same private mode-seeking primitives and are studied empirically in Section 5. Related work. Our work is closely related to a fast-growing literature on differentially private nonparametric methods. For density estimation, methods based on histograms, orthogonal series (Wasserman and Zhou, 2010), and kernels (Hall et al., 2013; Wagner et al., 2023; Liu et al., 2024) produce private density approximations. The closest prior work on private mode estimation is Pacchiano et al. (2021), who perturb a k-nearest-neighbor mode estimate and prove a differential privacy guarantee. In contrast, our procedure is a kernel-smoothed score-ascent method for recovering multiple density modes and is accompanied by smoothness-dependent upper and lower error rates. Turning to clustering, several works are based on the k-means framework (Balcan et al., 2017; Ghazi et al., 2020; Stemmer, 2021; Su et al., 2016), which optimize parametric objectives, but are not suited to nonparametric settings or irregular cluster shapes. Finally for regression tasks, most existing work on differential privacy focuses on approaches (Alabi et al., 2020; Arora et al., 2022; Cai et al., 2021, 2024; Wang, 2018; Sheffet, 2017) for modeling the conditional mean, leaving nonparametric modal regression unexplored. Finally, while differentially private mode estimation has been studied in Pacchiano et al. (2021) in the context of bandits, the authors consider a single mode and nearest neighbor based estimators are studied. Instead, we develop kernel based methods which lead to strictly improved rates when the density is sufficiently smooth. Note that nearest neighbor estimators are not equipped to take advantage of higher order smoothness. More importantly, we consider multiple modes, a feature that leads to significant complexity due to the inherent tension with privacy. While identifying each mode requires sufficiently granular data distribution, privacy requirements dictate that individual datum are still not distinguished. We therefore develop differentially private algorithms in a setting inspired by the study of mixture models: where the density modes are assumed to sufficiently separated and strongly identified. We also note the connection of our work with score-based denoising methods (see, for e.g., Ghosh et al., 2025; Wibisono et al., 2024). The need to preserve differential privacy in our case necessitates a combination between denoising step (score gradient ascent) and adding noise. While we use a specific kernel based score estimator, our approach can potentially be extended to the use of deep learning based score estimators and their differentially private versions. Organization. The remainder of the paper is organized as follows. Section 2 reviews some relevant background. Section 3 introduces the DP-GRAMS algorithm. Section 4 states the main theoretical results. Section 5 contains comprehensive empirical evaluation and implementation details. We conclude with a discussion of future work in Section 6. Proofs are deferred to Appendix A; additional experimental results are collected in Appendix B; and downstream pseudocode is given in Appendix C.

Background

This section reviews the ingredients underlying DP-GRAMS: population density modes, scorebased mode characterization, kernel density estimation, and differential privacy.

2.1

Modes as Statistical Objects

Our inferential target is the set of modes of a population density p : Rd → [0, ∞). A point µ ∈ Rd is a mode if it is a strict local maximizer of p. A standard sufficient second-order condition is ∇x p(x)|x=µ = 0

∇2x p(x)

and

x=µ

≺ 0.

For multimodal densities, the parameter of interest is therefore a finite, unordered collection of such critical points. We assume that p is strictly positive near each mode. Since the critical points of p coincide with those of ℓ(x) = log p(x), modes may equivalently be defined as local maximizers of log p. The gradient s(x) := ∇x log p(x) is the score function, and density modes are characterized by zeros of s(x) together with the corresponding local curvature condition. In the private algorithm, we work with a stabilized version of sb(x) defined below. All derivatives are taken with respect to the argument x; for notational simplicity, we omit the subscript x when no ambiguity arises.

2.2

Kernel Density Estimation

Given independent and identically distributed (i.i.d.) samples X1 , . . . , Xn ∈ Rd , the kernel density estimator (KDE, Chen (2017)) with bandwidth h > 0 is n 1 X x − Xi pb(x) = K , d nh i=1 h

where K is a kernel function integrating to one. The detailed conditions are stated in Section 4. Under standard regularity conditions, pb is smooth and admits well-defined derivatives. On regions where pb(x) > 0, it induces an estimator of the score function, ∇ log pb(x) =

∇pb(x) . pb(x)

To decrease the sensitivity of our score estimator from extreme observations, rather than working directly with ∇x log pb(x), we use a stabilized score estimator. Given a clipping level A > 0 and a density floor pfloor > 0, define clipA (∇pb(x)) sbA,pfloor (x) = , max{pb(x), pfloor }

A clipA (v) := v min 1, ∥v∥2

where

(2)

Whenever pb(x) ≥ pfloor and ∥∇pb(x)∥2 ≤ A, this coincides with the ordinary score estimator ∇ log pb(x). We focus on the behavior of pb and its derivatives near the population modes.

2.3

Mode Estimation via Gradient Ascent

To find critical points of the score function we employ a gradient ascent scheme on the score estimate: xt+1 = xt + η ∇ log pb(xt ). For Gaussian kernels, this gradient ascent on log pb matches exactly with the classical procedure of mean shift (Fukunaga and Hostetler, 1975; Cheng, 1995; Comaniciu and Meer, 2002). Unlike mean shift, the general gradient ascent formulation applies to arbitrary differentiable kernels. For our differentially private algorithm, we replace ∇ log pb(xt ) by the stabilized score estimator sbA,pfloor (xt ) and add Gaussian perturbations calibrated to its sensitivity. Differentiating the KDE, one has n 1 X x − Xi . ∇pb(x) = ∇K d+1 nh h i=1

Let us define

x − Xi . gi (x) := d+1 ∇K h h 1

Thus the stabilized score defined in (2) may be written as n 1 1X sbA,pfloor (x) = gi (x) . clipA max{pb(x), pfloor } n i=1

2.4

Differential Privacy

Differential privacy (Dwork et al., 2006, 2014) is a popular choice for guaranteeing protection of individual data. Definition 1 (Differential privacy (Dwork et al., 2014, Definition 2.4)). Let ε > 0 and δ ∈ [0, 1). A randomized mechanism M : X n → Y is (ε, δ)-differentially private if for all neighboring datasets X , X ′ differing in one entry and all measurable A ⊆ Y, Pr[M(X ) ∈ A | X ] ≤ eε Pr[M(X ′ ) ∈ A | X ′ ] + δ. We aim to find mode estimates that satisfy the above notion of privacy. To enforce such requirements, we use a stabilized score estimator together with Gaussian perturbations, placing our method within the framework of differentially private stochastic optimization. Our privacy analysis relies on standard tools including the Gaussian mechanism, privacy amplification by subsampling, composition, and post-processing invariance (see, e.g., Dwork et al., 2014; Balle et al., 2018; Dwork et al., 2010).

Differentially Private Mode Seeking Algorithm

In this section, we introduce our differentially private mode estimation algorithm. At a high level, the algorithm, called DP-GRAMS, proceeds in three stages: (i) it generates multiple initial points 6

concentrated in high-density regions with privacy guarantees, (ii) it refines each initialization via noisy ascent on the estimated score function, and (iii) it merges nearby outputs to produce a final set of private modes. Privacy and statistical guarantees for the full pipeline are established in Section 4. The complete procedure is summarized in Algorithms 1 and 2. Throughout, we write ε = εinit + εmodes , where εinit is the privacy budget allocated to initialization and εmodes to the gradient ascent stage. We now describe the different stages of this algorithm in more detail. Algorithm 1: DP-GRAMS: Differentially Private GRadient Ascent for Mode Seeking Input : Private data S = {Xi }ni=1 ⊂ Rd ; fixed public finite candidate set d cand Z = {zj }N j=1 ⊂ R ; privacy parameters (ε, δ) with ε = εinit + εmodes ; minibatch size m; number of iterations T ; number of DAP draws k ∈ N; DAP bandwidth hDAP > 0; suppression radius ρinit > 0; ascent bandwidth hmode > 0; stepsize η > 0; gradient clipping level A > 0; density floor pfloor > 0. c Output: Final private mode estimator M 1 Generate private initializations I = {x0,1 , . . . , x0,k } using Algorithm 2 with candidate set Z, k draws, bandwidth hDAP , and suppression radius ρinit ; write aℓ = x0,ℓ for the sampled anchors; 2 Compute the noise scale σ according to (9), using h = hmode ; 3 for t = 0, . . . , T − 1 do 4 Sample a minibatch Bt ⊂ [n] uniformly without replacement, with |Bt | = m; 5 Compute the current correlation matrix

(Kt )ℓr = C̄hmode (xt,ℓ , xt,r ) = exp − 6 7 8

1 ≤ ℓ, r ≤ k

Draw Zt ∈ Rk×d with independent columns (Zt )·,j ∼ N (0, σ 2 Kt ), for ℓ = 1, . . . , k do Compute the stabilized minibatch score sbA,pfloor ;Bt (xt,ℓ ) =

∥xt,ℓ − xt,r ∥ , hmode

j = 1, . . . , d.

clipA (∇pbBt (xt,ℓ )) max{pbBt (xt,ℓ ), pfloor }

Update

xt+1,ℓ = xt,ℓ + η sbA,pfloor ;Bt (xt,ℓ ) + (Zt )ℓ,·

f = {xT,1 , . . . , xT,k }; 10 Set M f to obtain M; c 11 Merge nearby points in M c 12 return M

Initialization. A central challenge in differentially private estimation of multiple modes is generating initial points that both respect privacy and lie in the modal basins of the true modes. To 7

Algorithm 2: Density-Aware Private (DAP) Initialization d cand Input : Private data S = {Xi }ni=1 ⊂ Rd ; public candidate set Z = {zj }N j=1 ⊂ R (see Section 4.2.2); privacy budget εinit ; number of DAP draws k ∈ N; DAP bandwidth hDAP > 0; suppression radius ρinit > 0. Output: Private initialization set I = {x0,1 , . . . , x0,k } 1 Set εdraw = εinit /k, A1 = [Ncand ]. 2 for j = 1, . . . , Ncand do 3 Compute n 1X uj = 1{∥Xi − zj ∥ ≤ hDAP }. n i=1 4 for ℓ = 1, . . . , k do

if Aℓ = ∅ then set Aℓ = [Ncand ];

Sample Jℓ ∈ Aℓ using

nεdraw Pr(Jℓ = j | Aℓ ) ∝ exp uj 1{j ∈ Aℓ }, 2

8 9

j ∈ [Ncand ].

Set aℓ = zJℓ and x0,ℓ = aℓ ; Update Aℓ+1 = Aℓ \ j ∈ [Ncand ] : ∥zj − aℓ ∥ ≤ ρinit .

10 return I = {x0,1 , . . . , x0,k };

address this, we use a Density-Aware Private (DAP) initialization scheme in Algorithm 2, which selects initializations from high-density regions while ensuring privacy via the exponential mechanism (McSherry and Talwar, 2007). To ensure that all high density neighborhoods are visited, we combine this algorithm with a local suppression step. cand More specifically, over a fixed public candidate set Z = {zj }N j=1 , we define the local empiricalmass utility n 1X uj = 1{∥Xi − zj ∥ ≤ hDAP }, j = 1, . . . , Ncand , (3) n i=1 which approximates the local probability mass near zj at scale hDAP . Since the sensitivity of uj is 1/n, the exponential mechanism uses weights proportional to nεdraw exp uj , 2

εdraw = εinit /k.

The anchor draws are performed one at a time. At round ℓ, Aℓ denotes the current active candidate set after suppression. The utility uj makes the scheme density-aware, since candidates with larger local empirical mass receive larger exponential-mechanism weight. After each selected anchor, we 8

remove candidate points in a ρinit -neighborhood of that anchor before the next draw, encouraging the algorithm to visit new high density regions and thus find previously unexplored modes. If the active set becomes empty before all k rounds are completed, the full candidate set is reopened. The selected candidate locations are used directly as the initialization set I = {x0,1 , . . . , x0,k }. The privacy guarantee for this stage is given by Theorem 4.1. In Section 4, we show that with a public hDAP -grid and k ≍ M log n draws, DAP places at least one initialization in each modal basin with high probability. Similar density-aware initialization ideas have been used earlier in Rodriguez and Laio (2014); Su et al. (2016); Li et al. (2016); Fan et al. (2023). In the rest of this paper, we write a1 , . . . , ak for the sampled anchors, so x0,ℓ = aℓ for each ℓ. Privatized Gradient Ascent. For the ascent stage, we use a possibly different bandwidth hmode , chosen according to (11). Given the KDE pb, DP-GRAMS iteratively updates candidate modes via ascent on an estimated score function. For any finite dataset X = {X1 , . . . , Xn } ⊂ Rd , define pbX (x) =

n X

nhdmode i=1

x − Xi , hmode

∇pbX (x) =

n X

nhd+1 mode i=1

∇K

x − Xi . hmode

(4)

Rather than working directly with the ordinary score ∇ log pbX (x), we use the stabilized score estimator in (2). In DP-GRAMS, the ascent update at round t is built from the minibatch field sbA,pfloor ;Bt (x) evaluated at the current iterates. Each candidate in the initialization pool I undergoes T iterations of minibatch ascent on the stabilized score function built using the bandwidth hmode . In each iteration, a minibatch Bt of size m is sampled uniformly without replacement, the minibatch score sbA,pfloor ;Bt (xt,ℓ ) is computed at each current iterate, and the resulting vector field is perturbed by Gaussian noise. Candidate modes are then updated using the fixed stepsize η > 0. After completing T iterations, the algorithm produces a set of private candidate modes f = {xT,1 , . . . , xT,k }. M Correlated Noise. The use of multiple initializations means that the same data are reused across several ascent trajectories. To account for this, at each iteration we use noise vectors that are independent across iterations and coordinates, but correlated across initializations. More specifically, at round t, after the current iterates xt,1 , . . . , xt,k are determined, we construct the correlation matrix ∥xt,ℓ − xt,r ∥ (Kt )ℓr = C̄hmode (xt,ℓ , xt,r ) = exp − , hmode

1 ≤ ℓ, r ≤ k,

based on the exponential kernel. The noise matrix Zt ∈ Rk×d then has independent columns with covariance σ 2 Kt . Thus trajectories whose current iterates are nearby receive more strongly correlated perturbations than trajectories whose current iterates are far apart. This is the correlated-noise mechanism used throughout DP-GRAMS. Merging Candidate Modes. Noise in gradient updates and randomness in initialization can produce multiple points corresponding to the same population mode. To reduce this redundancy, DP-GRAMS applies a final post-processing merge. When the number of modes is unknown, we use a radius-based merge with radius hmode , grouping candidate modes within distance hmode and 9

replacing each group by its mean. When the number of modes is known, we use Ward-linkage agglomerative clustering with target number of clusters equal to that mode count, again replacing c each final cluster by its mean. The final merged output is denoted by M.

3.1

Example: Modes of a Bivariate Gaussian u scores

DAP Initialization

6 Estimated KDE Surface for 4-Modal Bivariate Gaussian Mixture True modes

Density

0.020 0.015 0.010 0.005 0.000 5 5

0.125

0.100

0.075

0.050 0.025 5.0

Mode Estimation for 4-Modal Bivariate Gaussian Mixture True modes MS DP-GRAMS

0.150

x10

2.5

0.0

2.5

DP-GRAMS n=700 MS n=700 DP-GRAMS n=1000 MS n=1000 DP-GRAMS n=2000 MS n=2000 DP-GRAMS n=5000 MS n=5000

MSE

5.0

Privacy-Utility Tradeoff for 4-Modal Bivariate Gaussian Mixture

0.175

Data Candidate grid DAP Initializations

2 4 6 8

Privacy budget

Figure 1: (a) Estimated KDE surface with true modes overlaid. (b) DAP initialization: candidate points colored by the local empirical-mass utility uj in (3), together with the privately selected anchors used as the starting points for DP-GRAMS. (c) Contour plot comparing true (green), mean shift (blue), and DPGRAMS (orange) mode estimates on a single dataset. (d) Privacy–utility tradeoff: MSE vs. ε on a log scale for n ∈ {700, 1000, 2000, 5000} and ε ∈ {0.1, 0.25, 0.5, 1, 5}. Curves report averages over 20 runs with standard-error bars; dashed lines show non-private mean-shift baselines and solid curves show DP-GRAMS.

We use this example to visualize the initialization and final estimates from DP-GRAMS in a clean, well-separated modal landscape. Data are generated from the four-component bivariate Gaussian mixture (X, Y ) ∼ 14

4 X

N (µk , I2 ),

µ1 = (3, 3), µ2 = (3, −3), µ3 = (−3, 3), µ4 = (−3, −3),

k=1

which has four population modes at the corners of a square. Panels (a)–(c) are produced from one representative dataset with n = 1200 and (ε, δ) = (1, 10−6 ). The bandwidth is chosen by Silverman’s rule. 10

Figure 1 shows that the DAP initialization concentrates anchors in the four high-density basins, while the suppression step discourages redundant selections from the same local neighborhood. After private ascent, the final DP-GRAMS estimates remain close to both the population modes and the non-private mean-shift outputs. The privacy–utility curves show the largest gains when moving from very strict privacy to moderate privacy budgets, with further improvement as n increases. At the largest sample sizes and privacy budgets, the private estimates are close to the non-private mean-shift baseline in this well-separated benchmark.

Privacy Guarantees and Estimation Error Bounds

This section develops the privacy and utility guarantees for DP-GRAMS. We first state the regularity conditions under which the KDE and its derivatives are well behaved near the population modes. We then analyze privacy for the initialization and ascent stages, derive a basinwise error bound from a good initialization, and finally combine that local control with the DAP coverage argument to obtain global recovery of the final merged estimator. We first formalize the appropriate notions of local smoothness of the true density, and introduce higher-order kernels. For any r > 0, we write the Euclidean ball as B(µ, r) := {x ∈ Rd : ∥x−µ∥ ≤ r}. Definition 2 (Local Hölder class, (Definition 1.2 in Tsybakov, 2008)). Let U ⊂ Rd be open, and let β, L > 0. A function f : U → R belongs to Σ(β, L; U) if f is ⌊β⌋-times continuously differentiable on U and D⌊β⌋ f (x) − D⌊β⌋ f (x′ ) ≤ L ∥x − x′ ∥β−⌊β⌋ , x, x′ ∈ U. We say f is locally β-Hölder around µ on radius r if some open U satisfies B(µ, r) ⊂ U and f ∈ Σ(β, L; U ). Definition 3 (Kernel of order ℓ, (Definition 1.3 in Tsybakov, 2008)). A function K : Rd → R is a kernel of order ℓ ≥ 1 if Z

K(u) du = 1,

uα K(u) du = 0

for all 1 ≤ |α| ≤ ℓ,

and uα K(u) is integrable for all |α| ≤ ℓ. Before presenting our theoretical results, we present the required assumptions. We begin with regularity conditions on the smoothing kernel, which will be used to control KDE bias, stochastic fluctuation, and the argument underlying the correlated-noise mechanism. Assumption 1 (Kernel regularity). Let β > 2, ℓ = ⌊β⌋, and let K : Rd → R be a kernel of order ℓ as defined in Definition 3. Define K∞ := ∥K∥∞ and GK := supu∈Rd ∥∇K(u)∥. Assume: (i) supu∈Rd |∂ α K(u)| < ∞ for all multi-indices α with |α| ≤ 3. (ii) For all multi-indices α with |α| ≤ 2, (∂ α K)2 ∈ L1 (Rd ) and ∥u∥β |∂ α K(u)| ∈ L1 (Rd ); also, ∥∇K∥4 ∈ L1 (Rd ). 11

(iii) For all multi-indices α, α′ , γ with |α| ≤ 2, α′ ≤ α, and |γ| ≤ ℓ, uγ ∂ α K(u) ∈ L1 (Rd ) and ′ lim∥u∥→∞ uγ ∂ α K(u) = 0. (iv) ∂ α K ∈ L1 (Rd ) for all multi-indices α with |α| ≤ d + 2. We next impose local positivity and smoothness conditions on the population density in neighborhoods of the true modes. Assumption 2 (Model assumptions). Assume β > 2. Let p : Rd → [0, ∞) be a density taking finite values, with exactly M local modes at distinct points µ1 , . . . , µM such that min ∥µi − µj ∥ > c0 i̸=j

for a constant c0 > 0. Moreover, for each j ∈ [M ], there exist rj > 0, an open neighborhood Uj with B(µj , rj ) ⊂ Uj , and Lj > 0 such that: (i) pmin,j := inf x∈B(µj ,rj ) p(x) > c1 , pmax := supx∈Rd p(x) < c11 for a constant c1 > 0. (ii) p ∈ Σ(β, Lj ; Uj ), where Σ(β, Lj ; Uj ) is as defined in Definition 2. In particular, (ii) implies p ∈ C 2 (Uj ), and together with (i) the log-density ℓ(x) = log p(x) is well-defined on B(µj , rj ). Remark 4.1 (Spurious modes). Mode estimation is sensitive to outliers, which may create spurious local maxima in finite samples. Assumption 2(i) rules out population modes supported on vanishing mass by requiring p to be bounded away from zero near each µj . Algorithmically, even if an outlier is selected as an anchor, its influence is controlled by gradient clipping at level A together with the density floor pfloor . The final assumption in this preliminary block specifies the bandwidth regime used throughout the asymptotic analysis. Assumption 3 (Bandwidth condition). Let (hn )n≥1 be the bandwidth sequence. Assume hn ↓ 0

and

n hnd+4 →∞ log n

as n → ∞.

Remark 4.2 (Bandwidth notation). In the remainder of Section 4, unless stated otherwise, the bare bandwidth h refers to the ascent bandwidth hmode introduced in Section 3. The DAP initialization stage uses hDAP , which is written explicitly in the DAP construction, proposition, and proofs below. Together, these conditions yield uniform control of pb and its derivatives near modes, which we use for both privacy calibration and convergence.

4.1

Privacy and Sensitivity

We organize the privacy analysis in three steps. We first quantify the pointwise sensitivity of the stabilized score estimator, then calibrate the corresponding single-start Gaussian mechanism, and finally analyze the correlated multi-start mechanism used by DP-GRAMS. Lemma 1 (Deterministic sensitivity of the stabilized score). Assume Assumption 1. Let X = (X1 , . . . , Xn ) and X ′ = (X1′ , . . . , Xn′ ) be neighboring datasets, differing in one entry. Then, for every x ∈ Rd , Sh (A, pfloor ) , sbA,pfloor ;X (x) − sbA,pfloor ;X ′ (x) ≤ n where 2GK −(d+1) 2AK∞ −d h + 2 h . (5) Sh (A, pfloor ) := pfloor pfloor Lemma 1 identifies the relevant deterministic pointwise sensitivity of the stabilized score field, which directly yields a privacy calibration for a single ascent trajectory. Let us set δ δiter = , 2T

(

m εmodes εiter = min log 1 + (e − 1) , p , n 2 2T ln(2/δ)

εmodes 4T

)

For a single trajectory, we take Sh (A, pfloor )/m σ= log 1 + n(eεiter − 1)/m

2.5mT 2 log . nδ

(6)

The next lemma records the resulting privacy guarantee for a single privatized ascent run. Lemma 2 (Privacy of a single DP-GRAMS run). Consider Algorithm 1 with a fixed initialization x0 , stabilization parameters A and pfloor , and Gaussian noise scale σ given by (6). Then the final iterate xT is (εmodes , δ)-DP. Lemma 2 treats a single ascent trajectory started from one initialization. In DP-GRAMS, however, we run several trajectories in parallel from different anchors, and at each round these trajectories are evaluated on the same minibatch. Treating the k current score vectors as k separate Gaussian releases would lead to a sub-optimal privacy accounting, worse by a factor of the number of starts. To avoid this, we instead analyze the round-t collection of score vectors as a single joint Gaussian release whose covariance reflects the spatial proximity of the current iterates: nearby trajectories receive more strongly correlated perturbations, which allows the full multi-start release to be handled through a kernelized Gaussian mechanism in the spirit of Hall et al. (2013). Concretely, to couple the releases across starts without incurring an additional k-dependent privacy penalty, we use the exponential kernel ∥x − y∥ , C̄h (x, y) := exp − h b c for the Fourier transforms of C̄ := C̄1 and K, respectively. Define and write C̄ and K

I0 :=

1 (2π)d

Z Rd

c 2 |K(t)|

dt,

I1 :=

b C̄(t)

1 (2π)d

Z Rd

c 2 ∥t∥22 |K(t)| b C̄(t)

dt.

(7)

Under Assumption 1(iv), both I0 and I1 are finite. Let √ ∆h,corr (A, pfloor ) := 2 2

1/2

AI I1 h−(d+1) + 2 0 h−d . pfloor pfloor

(8)

The next lemma shows that the joint sensitivity of the concurrent stabilized score vectors is of the same order as the single-start sensitivity Sh (A, pfloor ), and therefore yields a joint privacy guarantee for the full multi-start ascent stage. Lemma 3 (Joint privacy across multiple initializations). Consider Algorithm 1 with a fixed initialization pool I = {x0,1 , . . . , x0,k }. For each round t = 0, . . . , T − 1, let us define Kt = [C̄h (xt,ℓ , xt,r )]kℓ,r=1

and

st (Bt ) := sbA,pfloor ;Bt (xt,ℓ )⊤ ℓ=1 ∈ Rkd .

Suppose that for each t, Yt = st (Bt ) + Ξt is the noisy joint score vector used to update the k trajectories, where Ξt ∼ N 0, σ 2 (Kt ⊗ Id ) are sampled independently across rounds, with ∆h,corr (A, pfloor )/m σ= log 1 + n(eεiter − 1)/m

2.5mT 2 log . nδ

(9)

Then the noisy score transcript (Y0 , . . . , YT −1 ) is (εmodes , δ)-DP. Consequently, the final merged c are (εmodes , δ)-DP. estimator M We next quantify the privacy cost of the data-dependent initialization pool produced by DAP. cand Theorem 4.1 (Privacy of DAP initialization). Let Z = {zj }N ⊂ Rd be a public grid fixed j=1 independently of the private sample X = (X1 , . . . , Xn ). Let the output of Algorithm 2 run for k rounds with total privacy budget εinit , be

I = {x0,1 , . . . , x0,k }. Then the initialization pool I is (εinit , 0)-DP. Combining the private initialization stage with the private correlated ascent stage yields the following end-to-end guarantee for DP-GRAMS. Corollary 4.2 (End-to-end privacy of DP-GRAMS). The complete DP-GRAMS algorithm, namely DAP Initialization with budget (εinit , 0), followed by the correlated multi-start ascent stage c is (εinit + εmodes , δ)-DP. of Lemma 3 with budget (εmodes , δ), and the final merge step producing M,

4.2

Error Bounds for Mode Estimation

We now turn to the utility analysis of DP-GRAMS. Since recovering multiple modes is a nonconvex problem (Carreira-Perpinan, 2007; Ota et al., 2019), the argument depends crucially on initialization. We therefore proceed in two steps. First, we study a single DP-GRAMS trajectory started in a local neighborhood of a true mode and derive its mean-squared error bound. Second, we show that the DAP design places at least one start in the local neighborhood of every population mode with high probability. 14

4.2.1

Local Convergence

Fix a mode µj for some j ∈ [M ]. We begin with the simpler problem in which one trajectory is already initialized inside the basin B(µj , rj ). Although the full algorithm uses correlated Gaussian perturbations across simultaneous starts, each individual row still has marginal law N (0, σ 2 Id ), so the basinwise analysis applies to any such trajectory. The next assumptions isolate the local geometry needed for this argument. Assumption 4 (Curvature at modes). Assume β > 2. For each j ∈ [M ] we assume that there exists a numerical constant C > 0 such that αj := −λmax ∇2 log p(µj ) > C.

Here, for a symmetric matrix A ∈ Rd×d , λmax (A) denotes its largest eigenvalue. We refer to αj as the local strong-concavity parameter of log p at µj . Under Assumption 2 with β > 2, ∇2 log p(x) is continuous, which together with Assumption 4 guarantees that, for each j ∈ [M ], there exists rej > 0 such that αj for all x ∈ B(µj , rej ). (10) ∇2 log p(x) ⪯ − Id 2 With this radius which guarantees sufficiently strong convexity of the Hessian, we define the local basins of attraction as follows. Definition 4 (Radius and separation). Assume β > 2. For each j ∈ [M ], let us define n

rj := min rej ,

o 1 min ∥µi − µj ∥ 2 i̸=j

where rej is as defined in (10). For each j ∈ [M ], we refer to B(µj , rj ) as the local basin neighborhood of µj . The above Definition 4 records the local-radius and separation conditions used in the following analysis. Note that by Assumption 2, the minimum separation between modes, i.e., mini̸=j ∥µi −µj ∥ is bounded below by a constant. Thus the definition of rj is not vacuous, and defines a strictly positive quantity. p Let us fix j ∈ [M ] and let cj := min,j 2 , where pmin,j is the local lower bound from Assumpb tion 2(i). Writing ℓ(x) = log p(x) and ℓ(x) = log pb(x) for convenience, we define a local good event as follows. The components of the local good event are: 1. The score function and its first two derivatives can be estimated suitably well by their sample counterparts. 2. The density estimate is bounded away from zero, while its variance and the squared norm of the kernel gradients remain bounded. 3. Once initialized in the basin of attraction of some mode, the privacy noise does not push the estimates out of that basin. 15

This event relates to the conditions that guarantee local convergence of our algorithm. Definition 5 (Local good event). Fix j ∈ [M ] and β > 2. Define the static local analytic event Alocal,j to be the event that there exist deterministic constants Cj , Dj > 0, depending only on the n local model and kernel constants and independent of n, such that 

sup

b − ∇s ℓ(x) ≤ Cj hβ−s + ∇s ℓ(x)

x∈B(µj ,rj )



log n  , nhd+2s

s = 0, 1, 2,

and n n X x − Xi 2 1 x − Xi 2 X ∇K , max K sup d h h x∈B(µj ,rj ) nh i=1 i=1 (

inf x∈B(µj ,rj )

pb(x) ≥ cj ,

)

≤ Dj .

Let us define the events stay,j := {xt ∈ B(µj , rj ) for all t = 0, . . . , T }, En,T

local,j stay,j := Alocal,j ∩ En,T . Gn,T n

and

local,j local,j Finally, let Xn,T be the set of all datasets X such that Gn,T holds.

The next theorem gives the basin-wise conditional MSE bound that drives the global recovery argument. The optimizing bandwidth is recorded immediately after the theorem. Theorem 4.3 (Local convergence of DP-GRAMS). Assume β > 2 and suppose Assumptions 1, 2, 3, 4 hold. Fix j ∈ [M ] and x0 ∈ B(µj , rj ). Then there exist constants pfloor , A, η̄j > 0 such that for every constant 0 < η ≤ η̄j , the following holds. Then, for the bandwidth h specified in (11), there local,j exist numerical constants CT , CnonDP,j , CDP,j > 0, such that for sufficiently large n, and Xn,T defined in Definition 5, local,j Pr X ∈ Xn,T ≥ 1 − 5n−4 , and i

E ∥xT − µj ∥2 X ≤ CnonDP,j

log n 2(β−1) d+2β

+ CDP,j

T d polylog(n, δ) β−1 d+β

n2 ε2modes

local,j . provided T = CT log n , m = ⌈n/ log n⌉ and X ∈ Xn,T

Remark 4.3. Some sufficient conditions for the choice of pfloor and A are given in Proposition A.1. The proof optimizes the underlying bandwidth-dependent error bound and yields the orderwise optimal choice  1 log n d+2β     n hopt ≍ 1  1 2d+2β    2 2

n εmodes

if εmodes ≳ εthr ,

εthr ≍ if εmodes ≲ εthr ,

log n n

−

d 2(d+2β)

1 n log n

1 2

(11)

When ε is large so the nonprivate term dominates, the resulting rate matches the minimax-optimal mode estimation rate of Arias-Castro et al. (2016); Genovese et al. (2014) up to logarithmic factors. We find that both the first and second terms of the MSE are affected by the traditional curse of dimensionality d, but is countered by the smoothness parameter β. If β → ∞, i.e., the density is analytic, we recover near-parametric rates. 16

4.2.2

Density-Aware Private Initialization

We now turn to DAP, whose role is to supply the good initializations required by the local convergence analysis. Recall that DAP is the initialization scheme in Algorithm 2, based on the utility in (3). We now specify the public candidate set used there. Let us consider a known public box Q=

d Y

[Lq , Uq ] ⊂ Rd ,

M [

such that

q=1

B(µj , rj ) ⊂ Q.

j=1

For each n, let us choose the number of DAP draws and the suppression radius so that ρinit ≍ (log n)−1/d .

k ≍ M log n,

Let Zn := Q ∩ hDAP Zd = {z1 , . . . , zNcand } be the candidate set used by Algorithm 2. We take

hDAP ≍

log n n

1 d+2β

(12)

d/(d+2β) , so that the public grid has polynomial size in n. Since Then Ncand ≍ h−d DAP ≍ (n/ log n) hDAP /ρinit → 0, the grid is fine enough at the suppression scale used below. The next proposition shows that, under the DAP design and a finite-sample lower bound on the initialization budget ϵinit , k ≍ M log n private draws localize in the modal basin neighborhoods and cover every such neighborhood with high probability.

Proposition 4.4 (High-probability coverage of DAP Initialization). Suppose Assumptions 2 and 4 hold, and Algorithm 2 is run under the DAP design above. Assume that εinit ≳ M n−2β/(d+2β) polylog(n)

and

k ≍ M log n.

Then there exists a constant Cinit,cov > 0, independent of εinit , such that, for all sufficiently large n,   PrI ⊆

M [

B(µj , rj ),

I ∩ B(µj , rj ) ̸= ∅ for every j ∈ [M ] ≥ 1 − Cinit,cov n−2 .

j=1

The above proposition pins together a number of crucial aspects for our initialization scheme. First of all, we ensure that the number of initializations scales as k ≍ M log n, which grows logarithmically in the sample size, and on average spends C log n initializations per mode. Secondly, on the high probability event defined within the proposition, the behavior of the chosen initializations is precisely characterized as follows. Each initialization lies within the neighborhood of some density mode; and, for each density mode, there exists at least one initialization in its convergence neighborhood. Both of these facts will be used in our global convergence analysis: since our iterative algorithm converges to local stationary points within the basin, which along with the first fact and the local concavity assumption, are precisely the modes. The second fact ensures that every mode has a nearby initialization and is thus guaranteed to be recovered. 17

4.2.3

Global Control

The next theorem upgrades the basinwise local control to the final merged estimator by combining DAP coverage with the deterministic post-processing conditions used by the merge rule. For each j ∈ [M ], let us define the basin index set Ij := {ℓ ∈ [k] : x0,ℓ ∈ B(µj , rj )}. which defines the set of particular initialization indices that are close to the j-th mode µj . Similar to our analysis for the local convergence, we now define a global good event, which in addition to the similar sub-events earlier defined for local convergence, also requires the initializations to be situated in the basins of attraction of each mode. n

Definition 6 (Global good event). Let us define the event of good initialization as Einit := ∀j ∈ o

[M ], I ∩ B(µj , rj ) ̸= ∅ and Aglobal := n

∩ Einit . For each j ∈ [M ] and ℓ ∈ Ij , let us stay,j define En,T (ℓ) := {xt,ℓ ∈ B(µj , rj ) for all t = 0, . . . , T }. Then we define the global good event as M local,j j=1 An

global Gn,T := Aglobal ∩ n

M \ \

stay,j En,T (ℓ).

j=1 ℓ∈Ij global global holds. be the set of all datasets X such that Gn,T Let Xn,T

Theorem 4.5 (Global convergence of DP-GRAMS). Assume β > 2 and suppose Assumptions 1, c be the output of Algorithm 1. Then there exist constants Cg , CnonDP,j , CP,j > 0 2, 3, 4 hold. Let M global such that for Xn,T defined in Definition 6, global ) ≥ 1 − Cglobal n−2 Pr(X ∈ Xn,T

c := {µ b1 , . . . , µ bM } satisfies for all sufficiently large n, and a permutation π : [M ] → [M ] such that M h

bj − µπ(j) ∥2 X ≤ CnonDP,j E ∥µ

log n 2(β−1) d+2β

+ CDP,j

T d polylog(n, δ) β−1 d+β

n2 ε2modes

global for every j ∈ [M ], whenever X ∈ Xn,T . global c corresponds Thus, on the global good event Gn,T , each point in the final merged estimator M to some population mode µπ(j) for an unknown permutation π : [M ] → [M ]. More importantly, the estimators satisfy the same conditional MSE rate as in Theorem 4.3.

4.3

Minimax Lower Bound

We now investigate the optimality of the mode estimation errors achieved by our algorithm. While score estimation under smoothness assumptions have been studied in previous literature, the problem of quantifying the loss in estimation accuracy due to privacy requirements, remains unexplored. To this end, we derive minimax lower bounds for mode estimation under the smoothness assumptions, while constraining the estimators to satisfy (ε, δ) privacy. 18

Theorem 4.6. Let β > 1 and let Pβ (L) denote the class of densities p : Rd → R+ satisfying Assumptions 2, 4 and that ∇2 log p(x) has bounded singular values for all x. Then there exists a constant C > 0 such that inf

b − x0 (p) sup Ep x

2(β−1)

− d+2β

≥ Cn

−

+ C(nε)

2(β−1) d+β

b x∈T (n,ε,δ) p∈Pβ (L)

provided δ = o(n−1 ). Here T (n, ε, δ) is the set of all possible estimators based on a sample of size n and satisfying (ε, δ) differential privacy. A comparison with Theorem 4.5 reveals that the error rates achieved by our DP-GRAMS algorithm are nearly minimax optimal, with the upper bound being worse than the lower bounds from Theorem 4.6 only up to the logarithmic terms log n and log δ. Our lower bound construction is based on first constructing a pointwise lower bound for differentially private score estimation, which to our knowledge is novel and might be of independent interest. More specifically, let the score function be sp (x) = ∇ log p(x). Then there exists a constant C > 0 such that inf

sup Ep sb − sp ∞

2(β−1)

− d+2β

≥ Cn

−

+ C(nε)

2(β−1) d+β

b s∈Tn,ε,δ (·) p∈Pβ (L)

provided δ = o(n−1 ). Given this result, the lower bound for mode estimation follows using the bounded singular values of the Hessian. The score estimation lower bound in turn depends on standard constructions for density perturbations, followed by the contraction of total variation under privacy constraints, from Karwa and Vadhan (2018).

Experiments

To complement the theoretical results, we evaluate DP-GRAMS and its downstream extensions on synthetic and real datasets. Unless stated otherwise, privacy is calibrated to (ε, δ), the ascent stage uses T = ⌈log n⌉ iterations and minibatch size m = ⌈n/ log n⌉, and reported summaries are averages over 20 independent runs with standard-error bars.

5.1

Implementation choices and tuning

We first summarize the implementation choices shared across all experiments. Kernel choice. The theory allows general differentiable kernels, but the experiments use the Gaussian kernel as the default. We validate this choice in Section 5.2 by comparing Gaussian DPGRAMS with an order-4 implementation. The Gaussian kernel is also convenient because classical mean shift is equivalent to a scaled ascent step on log pb.

Bandwidth selection. We use either the theoretically motivated h = ((log n)/n)1/(d+6) , or Silverman’s procedure (Silverman, 2018) to fix a single bandwidth h per dataset and keep it fixed across privacy budgets. Clipping. In the theoretical development, the ascent update is written using a clipped gradient and a density floor. In implementation, however, the key quantity is the size of the per-sample score contributions. Lemma 10 shows that, on the local good event, these contributions are bounded with high probability at the scale h−(d+1) . We therefore parameterize clipping through the effective threshold (13) C∗ = (clip multiplier) h−(d+1) . The dimensionless clip multiplier is treated as a dataset-level hyperparameter and selected once per dataset using a small pilot grid search. We then fix the resulting clip multiplier, and hence C∗ , for all privacy budgets and all reported runs. The density floor and gradient clipping level in the theoretical pseudocode are not separately tuned or reported in the experiments. Error metric. For all mode- and centroid-estimation tasks, accuracy is measured by the meanb1 , . . . , µ bb}: squared matching error between a reference set {µ1 , . . . , µk } and an estimated set {µ k MSE :=

min

b M: |M|=min{k,b max{k, k} k} (j,ℓ)∈M

bℓ ∥2 , ∥µj − µ

(14)

b with no repeated indices where the minimum is taken over all one-to-one matchings M ⊆ [k] × [k] b The denominator max{k, k} b keeps the error on the same scale and cardinality |M| = min{k, k}. even when the estimated number of modes differs from the reference number. The optimal matching is computed using the Hungarian algorithm (see, e.g., Kuhn, 1955; Munkres, 1957).

5.2

Kernel Choice Validation

We compare Gaussian DP-GRAMS with DP-GRAMS using a compactly supported order-4 kernel. The order-4 kernel is constructed from the orthonormal Legendre expansion (see Proposition 1.3 in Tsybakov, 2008) and extended to R2 by a product construction. We evaluate both methods on the four-corners Gaussian mixture and measure accuracy by the MSE formula in (14). To isolate the effect of the kernel, both methods use the same bandwidth h, computed from the order-β bandwidth rule in (11), for each (n, ε) configuration. Figure 2 shows that the Gaussian implementation is competitive with the order-4 alternative across the displayed privacy budgets and sample sizes. The order-4 kernel does not provide a systematic empirical advantage in this benchmark, while the Gaussian kernel avoids the finitesample complications associated with sign-indefinite KDE estimates. We therefore use the Gaussian kernel as the default in the remaining experiments.

MSE

Gaussian vs Order-4 DP-GRAMS on 4-Modal Bivariate Gaussian Mixture

103

= 0.2, = 10 6

= 0.5, = 10 6

= 1.0, = 10 6

= 2.0, = 10 6

2 × 103 Sample size n

5 × 103 Gaussian DP-GRAMS

103

2 × 103 Sample size n

5 × 103

Order-4 DP-GRAMS-

Figure 2: Kernel-choice validation on the four-corners Gaussian mixture. Gaussian DP-GRAMS is compared with order-4 DP-GRAMS. Privacy is calibrated with δ = 10−6 , and panels correspond to ε ∈ {0.2, 0.5, 1, 2}. Points show mean MSE (14) with standard-error bars.

5.3

Differentially Private Mode Estimation on Simulated Data

We next evaluate DP-GRAMS on simulated bivariate mixtures with known population modes. In each setting, we compare the population modes, non-private mean shift, and DP-GRAMS using the same bandwidth, iteration budget, and final merging rule. Privacy is calibrated with δ = 10−6 and ε ∈ {0.1, 0.25, 0.5, 1, 5}, and results are averaged over 20 independent runs with standard-error bars. The four-corners Gaussian mixture in Section 3.1 provides a clean, well-separated benchmark, while the five-component t-mixture below tests recovery under heavier tails and unequal local scales. Additional diagnostics for both synthetic mode-estimation benchmarks are reported in Appendix B.1, including repeated-run contour grids, sensitivity analyses for the clipping threshold C∗ , minibatch size m, and step size η, and full MSE and runtime tables across (n, ε). 5.3.1

Bivariate five-modal t-mixture

To examine robustness beyond the Gaussian setting in Section 3.1, we also consider the bivariate five-modal t-mixture (X, Y ) ∼

5 X

πk tνk (µk , σk2 I2 ),

πk = 0.2,

k=1

with centers µ1 = (0, 0), µ2 = (6, 0), µ3 = (−6, 0), µ4 = (0, 6), µ5 = (0, −6), degrees of freedom ν = (15, 6, 10, 8, 20), and scales σ = (0.1, 0.9, 1.3, 1.0, 0.4). This design is more challenging because its components have different tail behavior and local spread, so the resulting peaks are less homogeneous than in the Gaussian benchmark. 21

Mode Estimation for 5-Modal Bivariate t-Mixture

Estimated KDE Surface for 5-Modal Bivariate t-Mixture True modes

Privacy-Utility Tradeoff for 5-Modal Bivariate t-Mixture

True modes MS DP-GRAMS

10.0 7.5

x10

5 10

5.0 MSE

2.5 x2

Density

0.020 0.015 0.010 0.005 0.000

DP-GRAMS n=700 MS n=700 DP-GRAMS n=1000 MS n=1000 DP-GRAMS n=2000 MS n=2000 DP-GRAMS n=5000 MS n=5000

0.0 2.5 5.0 7.5 10

Privacy budget

Figure 3: Private mode estimation on the five-component t-mixture. Panels (a) and (b) use one representative dataset with n = 1200 and (ε, δ) = (1, 10−6 ). (a) Estimated KDE surface with true modes overlaid. (b) Contour plot comparing true (green), mean shift (blue), and DP-GRAMS (orange) mode estimates. (c) Privacy–utility tradeoff: MSE in (14) versus ε on a log scale for n ∈ {700, 1000, 2000, 5000} and ε ∈ {0.1, 0.25, 0.5, 1, 5}; dashed lines denote mean-shift baselines and solid curves show DP-GRAMS.

Figure 3 shows that DP-GRAMS continues to recover the modal structure in this harder setting. Panels (a) and (b) show peaks of unequal height and spread, yet the private estimates remain near the five modal locations on a representative dataset. Panel (c) shows that the privacy cost is largest for small samples and tight privacy budgets, while the gap to the mean-shift baseline narrows as n and ε increase. This behavior is consistent with the design: heavier tails and heterogeneous scales make the modal basins less uniform than in the four-corners Gaussian example.

5.4

Differentially Private Modal Regression on Simulated Data

Modal regression targets the conditional modes of Y | X = x, rather than the conditional mean. This distinction is important when the conditional response distribution is multimodal: meanbased smoothers can average across distinct subpopulations and fail to represent a typical response value (see Figure 4a). Building on the private ascent framework of DP-GRAMS, we obtain a differentially private analogue of partial mean shift (PMS, see Algorithm 1 of Chen et al. (2016a)), denoted DP-PMS; full pseudocode is given in Algorithm 4 in Section C. In the experiments below, we use the fixed-design version of this problem: the predictor locations X are treated as public, and privacy is enforced for the responses Y conditional on those predictors. This is appropriate for the simulated designs considered here, where X lies in a known public domain and the privacy-sensitive quantity is the response distribution. If the predictors themselves were private, an additional privacy mechanism would be needed for the predictor locations, e.g., by releasing the modal curve on a public evaluation grid or by privatizing the predictor-side binning. We evaluate PMS and DP-PMS on simulated regression problems with known oracle conditional modes. If yb(x) denotes an estimated modal response at predictor value x, and M(x) denotes the population conditional mode set, we use the pointwise loss

L yb(x), M(x) :=

min

m∈M(x)

yb(x) − m .

(15)

The reported regression error is the average of this loss over the predictor locations returned by the procedure. LOWESS is included only as a qualitative baseline targeting the conditional mean rather than the conditional modes. Privacy is calibrated to (ε, δ)-differential privacy with δ = 10−5 and ε ∈ {0.1, 0.2, 0.5, 1.0}. 5.4.1

Sinusoidal two-component mixture

To study performance under smoothly varying nonlinear modal structure, we consider a sinusoidal two-component mixture. We draw predictors X ∼ Uniform(0, 1) independently for each component and generate responses with Gaussian noise σ = 0.15 via Y1 = 1.5 + 0.5 sin(3πX) + N (0, σ 2 ),

Y2 = 0.5 sin(3πX) + N (0, σ 2 ),

so that the conditional density has two smooth modal curves. We consider sample sizes n ∈ {200, 600, 1200, 2400} and compare PMS, DP-PMS, and LOWESS on the same simulated design. Modal Regression for Sinusoidal 2-Component Mixture

1.5

1.0

1.0 Y

2.0

Privacy-Utility Tradeoff for Sinusoidal DP Modal Regression

Differentially Private Modal Regression for Sinusoidal Mixture

2.5

0.5

0.0

0.0 Data PMS LOWESS

0.5 1.0

0.0

0.2

0.5 0.4

0.6

0.8

1.0

PMS n=200 PMS n=600 PMS n=1200 PMS n=2400 DP-PMS n=200 DP-PMS n=600 DP-PMS n=1200 DP-PMS n=2400

MSE

2.5

Data DP-PMS ( =1)

0.0

0.2

0.4

0.6

0.8

1.0

Privacy budget

Figure 4: Private modal regression on sinusoidal two-component mixture data. Panels (a) and (b) use one representative dataset with n = 500; panel (b) uses (ε, δ) = (1, 10−5 ). (a) PMS captures the two conditional modes, whereas LOWESS averages across the mixture components. (b) DP-PMS recovers the same twobranch modal structure under privacy. (c) Privacy–utility tradeoff: oracle MSE in (15) versus ε on a log scale for n ∈ {200, 600, 1200, 2400}; points show averages over 20 runs with standard-error bars, and dashed lines show the non-private PMS baselines.

Figure 4 illustrates the qualitative distinction between mean and mode targets. PMS tracks the two modal branches, whereas LOWESS smooths across them. The private estimator preserves the two-branch structure in the representative run. The aggregate curves show the largest gains away from smallest n and tightest privacy budgets; for larger samples, further increases in ε produce smaller reductions in oracle MSE, so the remaining error is closer to the non-private PMS level. Additional modal-regression diagnostics are reported in Appendix B.2, including a complementary three-component piecewise-constant design, full privacy–utility and runtime tables for both regression benchmarks, and sensitivity analyses for clipping threshold C∗ and minibatch size m.

5.5

Differentially Private Modal Clustering

Private mode estimation also leads naturally to a private clustering procedure. In mode-based clustering, the modes of the density act as cluster representatives, and observations are assigned 23

to nearby modal centers. Building on DP-GRAMS, we obtain DP-GRAMS-C by first privately releasing a set of candidate modes, then merging them into a final collection of private centers, and finally assigning observations deterministically to the released centers. Full pseudocode is given in Algorithm 3 in Section C. c not to a separately released labelThe privacy guarantee applies to the released center set M, ing of the original private sample. In the experiments below, nearest-center assignments are used only for evaluation, for example when computing ARI and NMI against known labels. Thus DPGRAMS-C should be viewed as a private prototype-release procedure, analogous to differentially private k-means methods that release private cluster centers (Su et al., 2016): a curator releases private cluster representatives, and labels are then obtained by post-processing from those released representatives. As a private clustering baseline, we use DP-k-Means via the implementation of Su et al. (2016); Holohan et al. (2019). Density-threshold private clustering methods, such as private DBSCAN (Bozdemir et al., 2021; Qiu and Yi, 2025), are also relevant but target a different notion of cluster structure. Since DP-GRAMS-C releases centers and is evaluated with centroid-based metrics, DP-k-Means provides the closest standard private baseline. We report results at the same nominal ε, noting that DP-k-Means satisfies pure ε-DP whereas DP-GRAMS-C is calibrated under approximate (ε, δ)-DP with δ = o(n−1 ). Across all clustering experiments, we evaluate methods using adjusted Rand index (ARI), normalized mutual information (NMI), and centroid mean-squared error. The centroid error is computed using the MSE formula (14) between released centers and reference centroids. For highdimensional datasets, a full DAP lattice can become computationally prohibitive because the number of candidate grid points grows rapidly with dimension. We therefore run the image and geneexpression experiments in fixed PCA representations; MNIST uses public auxiliary candidates in PCA space, whereas Cancer RNA-Seq uses the default public DAP grid in reduced PCA space. 5.5.1

Simulated Gaussian blobs

We first evaluate DP-GRAMS-C on a simulated four-component Gaussian-blobs benchmark, where the true cluster identities and centroids are known. We generate two-dimensional datasets with k = 4 clusters using make blobs from scikit-learn (Pedregosa et al., 2011), with cluster standard deviation 1.2, and vary the sample size over n ∈ {700, 1000, 2000, 5000}. All clustering and centroid comparisons are performed in standardized feature space. For DP-GRAMS-C, privacy is calibrated with δ = 10−6 and ε ∈ {0.1, 0.2, 0.5, 1, 5}. We compare non-private mean-shift clustering, DP-GRAMS-C, standard k-means, and DP-k-Means; the target number of clusters is treated as known when applying the final agglomerative merge. Figures 5 and 6 show that DP-GRAMS-C recovers the underlying cluster geometry well in this clean benchmark. The released centers remain close to the true centroids and to the non-private mean-shift solution. Across privacy budgets and sample sizes, ARI and NMI increase while centroid MSE decreases as either n or ε increases. In this benchmark, DP-GRAMS-C is competitive with 24

Clustering Comparison on Blobs Dataset MS Clustering

DP-GRAMS-C

KMeans

1 0

x1 Cluster 0

Cluster 1

Cluster Assignments & Centroids Cluster 2

Cluster 3

True means

DP-KMeans

Estimated modes

Figure 5: Private clustering on a four-component blobs dataset (n = 1000). Each panel shows cluster assignments, true centroids, and estimated centroids for mean shift, DP-GRAMS-C, k-means, and DP-kMeans, with private methods run at ε = 1. Privacy-Utility (ARI) Tradeoff across n for Blobs

n = 700 0.7

0.6

NMI

0.5

n = 700

DP-GRAMS-C DP-KMeans

0.6

n = 1000 DP-GRAMS-C DP-KMeans

DP-GRAMS-C DP-KMeans

0.5

0.4

n = 2000

n = 5000

n = 2000

n = 5000

n = 2000

0.7

0.7 0.6

NMI

ARI

n = 1000

DP-GRAMS-C DP-KMeans

0.5

0.6

n = 5000 DP-GRAMS-C DP-KMeans

MSE

ARI

n = 1000 DP-GRAMS-C DP-KMeans

Privacy-Utility (MSE) Tradeoff across n for Blobs

MSE

n = 700 DP-GRAMS-C DP-KMeans

0.7

Privacy-Utility (NMI) Tradeoff across n for Blobs

DP-GRAMS-C DP-KMeans

0.5

0.4

DP-GRAMS-C DP-KMeans

modes

(a) ARI.

modes

DP-GRAMS-C DP-KMeans

modes

(b) NMI.

modes

0 10

modes

Figure 6: Privacy–utility tradeoff for private clustering on blobs across n ∈ {700, 1000, 2000, 5000} and ε ∈ {0.1, 0.2, 0.5, 1, 5}. Each panel shows a 2 × 2 grid over sample size, plotting ARI, NMI, or centroid MSE versus ε on a log scale for DP-GRAMS-C and DP-k-Means. Points show averages over 20 runs with standard-error bars.

DP-k-Means and is often stronger at moderate privacy budgets. Additional numerical summaries and sensitivity analyses for C∗ and minibatch size m are reported in Appendix B.3. 5.5.2

MNIST

We next study DP-GRAMS-C on MNIST. This experiment is intended to test the method at larger sample size in a moderate-dimensional representation. Because MNIST is public, we construct a stratified public auxiliary candidate set of 1000 images and project these candidates into the same PCA space used for clustering. These points are used as DAP candidates, while the remaining MNIST images form the private experimental sample; no full DAP lattice grid is built for MNIST. We standardize the pixel features, fit a whitened five-dimensional PCA representation, and run clustering in this reduced space. For centroid-error evaluation, estimated centers are mapped back to standardized pixel space using the inverse PCA map and compared against standardized 25

class means computed on the private experimental sample. We compare four methods: non-private mean-shift clustering, DP-GRAMS-C, standard k-means, and DP-k-Means. For DP-GRAMS-C, privacy is calibrated with δ = 10−5 and ε ∈ {0.05, 0.1, 0.2, 0.5, 1}. Clustering Comparsion on MNIST Dataset DP-GRAMS-C

KMeans

DP-KMeans

0 2

PC1

PC2

MS Clustering

0 2

PC1

Cluster 0 Cluster 1

Cluster 2 Cluster 3

4 Cluster 4 Cluster 5

2 2

Cluster Assignments & Centroids Cluster 6 Cluster 7

Cluster 8 Cluster 9

PC1

True means Estimated modes

Figure 7: MNIST dataset with public auxiliary candidates. Two-dimensional visualization of the fivedimensional PCA clustering representation, comparing mean shift, DP-GRAMS-C, k-means, and DP-kMeans, with private methods run at ε = 1. True class centroids and estimated centroids are overlaid.

Privacy-Utility: ARI (MNIST)

Privacy-Utility: NMI (MNIST)

Privacy-Utility: MSE (MNIST)

0.425

0.30

0.400

300

NMI

ARI

MSE

0.375

0.25

0.350 0.325

0.20

0.300

DP-GRAMS-C DP-KMeans

0.15 10

modes

250 200 150

DP-GRAMS-C DP-KMeans

0.275 0

DP-GRAMS-C DP-KMeans

350

100 0

modes

Figure 8: Privacy–utility on MNIST with public auxiliary candidates in the whitened five-dimensional PCA representation: ARI, NMI, and centroid MSE versus ε on a log scale for DP-GRAMS-C and DP-k-Means, with ε ∈ {0.05, 0.1, 0.2, 0.5, 1}. Points show averages over 20 runs with standard-error bars.

Figures 7 and 8 show that DP-GRAMS-C performs well in the whitened five-dimensional MNIST PCA representation. The public auxiliary candidates avoid the combinatorial growth of a full DAP grid, while the released private centers retain class structure visible in the PCA visualization. The largest improvement occurs from ε = 0.05 to moderate privacy budgets; at ε = 1, DP-GRAMS-C is close to non-private mean shift in ARI and NMI and remains stronger than the DP-k-Means baseline. Additional MNIST summaries and sensitivity analyses for C∗ and minibatch size m are reported in Appendix B.4. 5.5.3

Cancer RNA-Seq

We also evaluate DP-GRAMS-C on the UCI Gene Expression Cancer RNA-Seq dataset of Fiorini (2016), which consists of high-dimensional RNA-seq gene-expression profiles with five tumor types (BRCA, COAD, KIRC, LUAD, PRAD). We standardize each gene to zero mean and unit variance, 26

project the data to a whitened six-dimensional PCA representation, and run all clustering methods in this reduced space. Unlike MNIST, this experiment uses the default public DAP grid in the reduced PCA space rather than public auxiliary candidate images. For centroid-error evaluation, estimated centers are mapped back to standardized gene space through the inverse PCA map and compared against standardized class means. We compare mean-shift clustering, DP-GRAMS-C, standard k-means, and DP-k-Means. For DP-GRAMS-C, privacy is calibrated with δ = 10−5 and ε ∈ {0.5, 1, 2, 5, 10}. Clustering Comparison on Cancer RNA-Seq Data

PC2

0 1

KMeans

PC2

DP-GRAMS-C

DP-KMeans

2 1 0

PC2

MS Clustering

PC1

Cluster 0 Cluster 1

Cluster 2

2 Cluster Assignments & Centroids Cluster 3

4 4

Cluster 4

True means

PC1

Estimated modes

Figure 9: Cancer RNA-Seq dataset after gene-wise standardization and projection to six principal components. Two-dimensional PCA visualization of the six-dimensional clustering representation, comparing mean shift, DP-GRAMS-C, k-means, and DP-k-Means, with private methods run at ε = 1. True class centroids and estimated centroids are overlaid.

Privacy-Utility: ARI (Cancer RNA-Seq)

Privacy-Utility: NMI (Cancer RNA-Seq)

0.8

0.4

0.6

DP-GRAMS-C DP-KMeans

MSE

NMI

0.6

ARI

Privacy-Utility: MSE (Cancer RNA-Seq)

0.4 DP-GRAMS-C DP-KMeans

0.2 10

DP-GRAMS-C DP-KMeans

0.2 10

Figure 10: Privacy–utility on Cancer RNA-Seq after gene-wise standardization and projection to six principal components: ARI, NMI, and centroid MSE versus ε on a log scale for DP-GRAMS-C and DP-k-Means, with ε ∈ {0.5, 1, 2, 5, 10}. Points show averages over 20 runs with standard-error bars.

Figures 9 and 10 show that the private mode-based procedure remains effective on the Cancer RNA-Seq task after reduction to a whitened six-dimensional PCA representation. The released centers preserve the class separation visible in the PCA display, and the aggregate curves show a sharp improvement from ε = 0.5 to ε = 1: ARI and NMI move near the non-private range, while centroid MSE decreases as the privacy budget increases. Table 1 provides a cross-dataset one-run comparison at ε = 1, where DP-GRAMS-C is close to non-private mean shift in ARI and NMI and much stronger than DP-k-Means on Cancer RNA-Seq. Additional Cancer RNA-Seq 27

Table 1: One-run clustering performance at nominal privacy budget ε = 1 across datasets. MSE denotes centroid mean-squared error. Blobs NMI

MNIST MSE

ARI

NMI

Cancer RNA-Seq

Algorithm

ARI

MSE

ARI

NMI

MSE

Mean shift

0.760 0.720 0.0437 0.314 0.424 113.294 0.816 0.846

527.335

DP-GRAMS-C 0.727 0.704 0.0562 0.331 0.433

53.489

0.797 0.845

700.850

k-Means

0.753 0.715 0.0018 0.295 0.409

68.104

0.817 0.863

340.236

DP-k-Means

0.620 0.638 0.2904 0.302 0.411 176.529 0.054 0.120 14822.003

summaries, including full privacy–utility tables and sensitivity analyses for C∗ and minibatch size m, are reported in Appendix B.5.

Discussion

This paper introduced DP-GRAMS, a differentially private mode-seeking algorithm that leverages the equivalence between mean shift and gradient ascent on the log-density, bringing differentially private stochastic optimization tools to nonparametric mode estimation. The key insight is to decompose the log-KDE score field into per-sample contributions, clip those contributions to control sensitivity, and add Gaussian noise calibrated via standard (ε, δ)-DP accounting. This perspective enables principled privatization of a classical nonparametric algorithm while preserving its underlying geometric interpretation. Our analysis points to multiple directions that might be of interest. Firstly, initialization plays a crucial role in a mode seeking problem, especially in settings where local modal basins can be separated by low-density regions. To address this under privacy constraints, we incorporate a simple differentially private initialization scheme that combines a density-aware exponentialmechanism utility with a local suppression step. In our analysis, the public hDAP -grid together with the suppression rule yields a private initialization scheme that visits every modal basin with high probability using only k ≍ M log n draws. Further improving the efficiency of private initialization might be a promising direction. Secondly, from a privacy accounting perspective, we adopt the standard (ε, δ)-differential privacy framework rather than alternatives such as Rényi DP Mironov (2017) or zero-Concentrated DP Bun and Steinke (2016). While RDP- and zCDP-based analyses can yield tighter composition bounds in some regimes, (ε, δ)-DP remains a widely used notion with a direct and interpretable guarantee. Investigating RDP or zCDP variants for score and mode estimationmay be another avenue for future work. Finally, there is potential in adapting our results to unknown smoothness and additional structures that the true density might enjoy. Adapting to the smoothness β, with and without privacy is interesting: see, e.g., Lepskii (1991); Kroll (2019); Butucea et al. (2020); Schluttenhofer

and Johannes (2022); Auddy et al. (2025). More interestingly, deep learning based estimators have recently shown immense promise in adapting to underlying dimensionality and specific dependence patterns of the score function: see, e.g., Nakada and Imaizumi (2020); Song and Ermon (2019); Oko et al. (2023). Advancing these results to incorporate differential privacy is of both theoretical and practical interest. We intend to pursue this in the future.

References Daniel Alabi, Audra McMillan, Jayshree Sarathy, Adam Smith, and Salil Vadhan. Differentially private simple linear regression. arXiv preprint arXiv:2007.05157, 2020. Ery Arias-Castro, David Mason, and Bruno Pelletier. On the estimation of the gradient lines of a density and the consistency of the mean-shift algorithm. The Journal of Machine Learning Research, 17(1):1487–1514, 2016. Raman Arora, Raef Bassily, Cristóbal Guzmán, Michael Menart, and Enayat Ullah. Differentially private generalized linear models revisited. Advances in neural information processing systems, 35:22505–22517, 2022. Arnab Auddy, T Tony Cai, and Abhinav Chakraborty. Minimax and adaptive transfer learning for nonparametric classification under distributed differential privacy constraints. Journal of the Royal Statistical Society Series B: Statistical Methodology, page qkaf070, 2025. Shai Avidan. Ensemble tracking. IEEE transactions on pattern analysis and machine intelligence, 29(2):261–271, 2007. Maria-Florina Balcan, Travis Dick, Yingyu Liang, Wenlong Mou, and Hongyang Zhang. Differentially private clustering in high-dimensional euclidean spaces. In International Conference on Machine Learning, pages 322–331. PMLR, 2017. Borja Balle, Gilles Barthe, and Marco Gaboardi. Privacy amplification by subsampling: Tight analyses via couplings and divergences. Advances in neural information processing systems, 31, 2018. Raef Bassily, Adam Smith, and Abhradeep Thakurta. Differentially private empirical risk minimization: Efficient algorithms and tight error bounds. arXiv preprint arXiv:1405.7085, 2014. Beyza Bozdemir, Sébastien Canard, Orhan Ermis, Helen Möllering, Melek Önen, and Thomas Schneider. Privacy-preserving density-based clustering. In Proceedings of the 2021 ACM Asia Conference on Computer and Communications Security, pages 658–671, 2021. Mark Bun and Thomas Steinke. Concentrated differential privacy: Simplifications, extensions, and lower bounds. In Theory of cryptography conference, pages 635–658. Springer, 2016.

Cristina Butucea, Amandine Dubois, Martin Kroll, and Adrien Saumard. Local differential privacy: Elbow effect in optimal density estimation and adaptation over besov ellipsoids. Bernoulli, 26 (3):1727–1764, 2020. T Tony Cai, Yichen Wang, and Linjun Zhang. The cost of privacy: Optimal rates of convergence for parameter estimation with differential privacy. The Annals of Statistics, 49(5):2825–2850, 2021. T Tony Cai, Abhinav Chakraborty, and Lasse Vuursteen. Optimal federated learning for nonparametric regression with heterogeneous distributed differential privacy constraints. arXiv preprint arXiv:2406.06755, 2024. Miguel A Carreira-Perpinan. Gaussian mean-shift is an em algorithm. IEEE Transactions on Pattern Analysis and Machine Intelligence, 29(5):767–776, 2007. José E Chacón. A population background for nonparametric density-based clustering. Statistical Science, 30(4):518–532, 2015. Yen-Chi Chen. A tutorial on kernel density estimation and recent advances. Biostatistics & Epidemiology, 1(1):161–187, 2017. Yen-Chi Chen, Christopher R. Genovese, Ryan J. Tibshirani, and Larry Wasserman. Nonparametric modal regression. The Annals of Statistics, 44(2):489 – 514, 2016a. doi: 10.1214/15-AOS1373. URL https://doi.org/10.1214/15-AOS1373. Yen-Chi Chen, Christopher R. Genovese, and Larry Wasserman. A comprehensive approach to mode clustering. Electronic Journal of Statistics, 10(1):210–241, 2016b. doi: 10.1214/15-EJS1102. Yizong Cheng. Mean shift, mode seeking, and clustering. IEEE transactions on pattern analysis and machine intelligence, 17(8):790–799, 1995. Dorin Comaniciu and Peter Meer. Mean shift: A robust approach toward feature space analysis. IEEE Transactions on pattern analysis and machine intelligence, 24(5):603–619, 2002. Dorin Comaniciu, Visvanathan Ramesh, and Peter Meer. Kernel-based object tracking. IEEE Transactions on pattern analysis and machine intelligence, 25(5):564–577, 2003. Fida Kamal Dankar and Khaled El Emam. Practicing differential privacy in health care: A review. Trans. Data Priv., 6(1):35–67, 2013. Travis Dick, Cynthia Dwork, Michael Kearns, Terrance Liu, Aaron Roth, Giuseppe Vietri, and Zhiwei Steven Wu. Confidence-ranked reconstruction of census microdata from published statistics. Proceedings of the National Academy of Sciences, 120(8):e2218605120, 2023. Cynthia Dwork. Differential privacy. In International colloquium on automata, languages, and programming, pages 1–12. Springer, 2006. 30

Cynthia Dwork, Krishnaram Kenthapadi, Frank McSherry, Ilya Mironov, and Moni Naor. Our data, ourselves: Privacy via distributed noise generation. In Annual international conference on the theory and applications of cryptographic techniques, pages 486–503. Springer, 2006. Cynthia Dwork, Guy N Rothblum, and Salil Vadhan. Boosting and differential privacy. In 2010 IEEE 51st annual symposium on foundations of computer science, pages 51–60. IEEE, 2010. Cynthia Dwork, Aaron Roth, et al. The algorithmic foundations of differential privacy. Foundations and Trends® in Theoretical Computer Science, 9(3–4):211–407, 2014. Jochen Einbeck and Gerhard Tutz. The fitting of multifunctions: an approach to nonparametric multimodal regression. A. Rizzi, MV, editor, COMPSTAT 2006, Proceedings in Computational Statistics, pages 1243–1250, 2006. Chenglin Fan, Ping Li, and Xiaoyun Li. K-median clustering via metric embedding: towards better initialization with differential privacy. Advances in Neural Information Processing Systems, 36: 73817–73838, 2023. Samuele Fiorini. gene expression cancer RNA-Seq. UCI Machine Learning Repository, 2016. DOI: https://doi.org/10.24432/C5R88H. Keinosuke Fukunaga and Larry Hostetler. The estimation of the gradient of a density function, with applications in pattern recognition. IEEE Transactions on information theory, 21(1):32–40, 1975. Christopher R. Genovese, Marco Perone-Pacifico, Isabella Verdinelli, and Larry Wasserman. Nonparametric ridge estimation. The Annals of Statistics, 42(4):1511 – 1545, 2014. doi: 10.1214/14-AOS1218. URL https://doi.org/10.1214/14-AOS1218. Badih Ghazi, Ravi Kumar, and Pasin Manurangsi. Differentially private clustering: Tight approximation ratios. Advances in Neural Information Processing Systems, 33:4040–4054, 2020. Sulagna Ghosh, Nikolaos Ignatiadis, Frederic Koehler, and Amber Lee. Stein’s unbiased risk estimate and hyv\” arinen’s score matching. arXiv preprint arXiv:2502.20123, 2025. Rob Hall, Alessandro Rinaldo, and Larry Wasserman. Differential privacy for functions and functional data. The Journal of Machine Learning Research, 14(1):703–727, 2013. Naoise Holohan, Stefano Braghin, Pól Mac Aonghusa, and Killian Levacher. Diffprivlib: the ibm differential privacy library. arXiv preprint arXiv:1907.02444, 2019. Vishesh Karwa and Salil Vadhan. Finite sample differentially private confidence intervals. In 9th Innovations in Theoretical Computer Science Conference (ITCS 2018), volume 94, page 44. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, 2018.

Martin Kroll. Pointwise adaptive kernel density estimation under local approximate differential privacy. arXiv preprint arXiv:1907.06233, 2019. Harold W Kuhn. The hungarian method for the assignment problem. Naval research logistics quarterly, 2(1-2):83–97, 1955. OV Lepskii. On a problem of adaptive estimation in gaussian white noise. Theory of Probability & Its Applications, 35(3):454–466, 1991. Fengfu Li, Hong Qiao, and Bo Zhang. Effective deterministic initialization for k-means-like methods via local density peaks searching. arXiv preprint arXiv:1611.06777, 2016. Jia Li, Surajit Ray, and Bruce G Lindsay. A nonparametric statistical approach to clustering via mode identification. Journal of Machine Learning Research, 8(Aug):1687–1723, 2007. Erzhi Liu, Jerry Yao-Chieh Hu, Alex Reneau, Zhao Song, and Han Liu. Differentially private kernel density estimation. arXiv preprint arXiv:2409.01688, 2024. 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. Ilya Mironov. Rényi differential privacy. In 2017 IEEE 30th computer security foundations symposium (CSF), pages 263–275. IEEE, 2017. James Munkres. Algorithms for the assignment and transportation problems. Journal of the society for industrial and applied mathematics, 5(1):32–38, 1957. Ryumei Nakada and Masaaki Imaizumi. Adaptive approximation and generalization of deep neural network with intrinsic dimensionality. Journal of Machine Learning Research, 21(174):1–38, 2020. Daniel L Oberski and Frauke Kreuter. Differential privacy and social science: An urgent puzzle. Harvard Data Science Review, 2(1):1–21, 2020. Kazusato Oko, Shunta Akiyama, and Taiji Suzuki. Diffusion models are minimax optimal distribution estimators. In International Conference on Machine Learning, pages 26517–26582. PMLR, 2023. Hirofumi Ota, Kengo Kato, and Satoshi Hara. Quantile regression approach to conditional mode estimation. Electronic Journal of Statistics, 13:3120–3160, 2019. Aldo Pacchiano, Heinrich Jiang, and Michael I Jordan. Robustness guarantees for mode estimation with an application to bandits. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 9277–9284, 2021.

F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825–2830, 2011. Yuan Qiu and Ke Yi. Approximate dbscan under differential privacy. Proceedings of the ACM on Management of Data, 3(3):1–24, 2025. Alex Rodriguez and Alessandro Laio. Clustering by fast search and find of density peaks. science, 344(6191):1492–1496, 2014. Sandra Schluttenhofer and Jan Johannes. Adaptive pointwise density estimation under local differential privacy. arXiv preprint arXiv:2206.07663, 2022. Or Sheffet. Differentially private ordinary least squares. In International Conference on Machine Learning, pages 3105–3114. PMLR, 2017. Bernard W Silverman. Density estimation for statistics and data analysis. Routledge, 2018. Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems, 32, 2019. Uri Stemmer. Locally private k-means clustering. Journal of Machine Learning Research, 22(176): 1–30, 2021. Dong Su, Jianneng Cao, Ninghui Li, Elisa Bertino, and Hongxia Jin. Differentially private k-means clustering. In Proceedings of the sixth ACM conference on data and application security and privacy, pages 26–37, 2016. Aleksandr Borisovich Tsybakov. Recursive estimation of the mode of a multivariate distribution. Problemy Peredachi Informatsii, 26(1):38–45, 1990. Alexandre B Tsybakov. Nonparametric estimators. In Introduction to Nonparametric Estimation, pages 1–76. Springer, 2008. Tal Wagner, Yonatan Naamad, and Nina Mishra. Fast private kernel density estimation via locality sensitive quantization. In International Conference on Machine Learning, pages 35339–35367. PMLR, 2023. Yu-Xiang Wang. Revisiting differentially private linear regression: optimal and adaptive prediction & estimation in unbounded domain. arXiv preprint arXiv:1803.02596, 2018. Larry Wasserman and Shuheng Zhou. A statistical framework for differential privacy. Journal of the American Statistical Association, 105(489):375–389, 2010.

Andre Wibisono, Yihong Wu, and Kaylee Yingxi Yang. Optimal score estimation via empirical bayes smoothing. In The Thirty Seventh Annual Conference on Learning Theory, pages 4958– 4991. PMLR, 2024. Bin Yu. Assouad, fano, and le cam. In Festschrift for Lucien Le Cam: research papers in probability and statistics, pages 423–435. Springer, 1997.

Proofs

This appendix proves the theoretical results in the order in which they are used. We begin with the privacy statements, then assemble the local analytic and probabilistic ingredients behind Theorem 4.3. We next prove DAP coverage and combine it with the local theorem and deterministic post-processing to prove Theorem 4.5. We then prove the minimax lower bound. The final subsection collects the remaining auxiliary proofs in the same order of invocation.

A.1

Privacy results

We first verify privacy one stage at a time. The initialization pool is handled through repeated exponential-mechanism draws with suppression, after which the end-to-end guarantee follows by composing that initialization step with the correlated Gaussian ascent mechanism. Proof of Theorem 4.1. Fix neighboring datasets X = (X1 , . . . , Xn ) and X ′ = (X1′ , . . . , Xn′ ) differing in one entry. Since the candidate set Z is fixed independently of the data, each score uj (X ) =

n 1X 1{∥Xi − zj ∥ ≤ hDAP } n i=1

has global sensitivity at most 1/n. Set εdraw = εinit /k. Fix a round ℓ ∈ [k], and condition on the previously selected indices J1 , . . . , Jℓ−1 . Under this conditioning, the set Aℓ is fixed. Hence the range from which Jℓ is sampled is also fixed. More precisely, define

Rℓ :=

   Aℓ ,

if Aℓ ̸= ∅,

 [Ncand ],

if Aℓ = ∅.

Then Algorithm 2 samples Jℓ from this fixed range Rℓ . Thus, for j ∈ [Ncand ], exp nεdraw uj (X ) 1{j ∈ Rℓ } . Pr(Jℓ = j | Rℓ , X ) = P 2 nεdraw r∈Rℓ exp 2 ur (X )

This is exactly the exponential mechanism on the fixed restricted range Rℓ with score sensitivity 1/n (McSherry and Talwar, 2007). Hence round ℓ is (εdraw , 0)-DP.

By sequential composition, the full index sequence (J1 , . . . , Jk ) is k X

εdraw , 0 = (εinit , 0)-DP

ℓ=1

(see, e.g., Dwork et al., 2014). Finally, since aℓ = zJℓ and x0,ℓ = aℓ , both the anchor sequence (a1 , . . . , ak ) and the initialization pool I = {x0,1 , . . . , x0,k } are deterministic post-processing of (J1 , . . . , Jk ), and are therefore also (εinit , 0)-DP (Dwork et al., 2014, Proposition 2.1). Proof of Corollary 4.2. By Theorem 4.1, the DAP initialization pool I = {x0,1 , . . . , x0,k } is (εinit , 0)-DP. Conditional on any fixed initialization pool I, Lemma 3 shows that the remainder c is (εmodes , δ)-DP. Therefore, by sequential of Algorithm 1, including the final merged estimator M, composition, (εinit , 0) + (εmodes , δ) = (εinit + εmodes , δ), so the complete DP-GRAMS algorithm is (εinit + εmodes , δ)-DP (see, e.g., Dwork et al., 2014).

A.2

Local convergence: ingredient statements and proof roadmap

Fix j ∈ [M ] and consider a single DP-GRAMS trajectory x0 , . . . , xT initialized in B(µj , rj ). In this local block, xt denotes this single trajectory and zt ∼ N (0, σ 2 Id ) denotes its Gaussian perturbation at round t. The argument starts from the local geometry of log p near µj , upgrades this to uniform control of the KDE and log-KDE on the modal basin, packages those bounds into the event Alocal,j , n and then combines the stopped recursion with basin retention before optimizing in h. Lemma 4 (Continuity of ∇2 log p at each mode). Under Assumption 2 with β > 2, for each j ∈ [M ], the map x 7→ ∇2 log p(x) is continuous at µj . Consequently, for every ξ > 0 there exists rej (ξ) > 0 such that ∇2 log p(x) − ∇2 log p(µj ) ≤ ξ,

∀x ∈ B(µj , rej (ξ)).

In particular, under Assumption 4, taking ξ = αj /2 yields ∇2 log p(x) ⪯ −

αj Id , 2

∀x ∈ B(µj , rej (αj /2)).

Lemma 5 (Derived inward drift on the local neighborhood). Assume β > 2 and suppose Assumptions 2, 4, and 4 hold. Then, for each j ∈ [M ] and every x ∈ B(µj , rj ), ⟨x − µj , ∇ log p(x)⟩ ≤ −

αj ∥x − µj ∥2 . 2

The next elementary geometric observation records the role of the radius restriction in Assumption 4: the chosen local basin neighborhood around µj contains no other population mode. 35

Lemma 6 (Uniform KDE derivative rates on local balls). Assume β > 2 and let ℓ = ⌊β⌋ ≥ 2. Let pb be the KDE defined in Section 3, based on i.i.d. data from a density p satisfying Assumptions 2 and 3, with kernel K satisfying Assumption 1. Fix j ∈ [M ]. Then with probability at least 1 − n−4 , simultaneously for all s ∈ {0, 1, 2}, 

sup

∥∇s pb(x) − ∇s p(x)∥ ≤ C hβ−s +

x∈B(µj ,rj )



log n  , nhd+2s

where C > 0 depends only on d, β, the local Hölder constants of p on B(µj , rj ), and kernel moments. Lemma 7 (Uniform lower bound for pb on mode basins). Let pb be the KDE from Section 3, built from i.i.d. samples drawn from p satisfying Assumptions 2, 3, and 1. Fix j ∈ [M ]. Then there exists a constant cj > 0 such that for all sufficiently large n, !

inf x∈B(µj ,rj )

pb(x) ≥ cj

≥ 1 − n−4 .

Lemma 8 (Uniform-in-x log-KDE derivative rates on local balls). Fix j ∈ [M ] and assume β > 2. Let ℓ = log p and ℓb = log pb. There exists a constant Cj > 0 (depending only on pmin,j and supx∈B(µj ,rj ) ∥∇r p(x)∥ for r ≤ 2) such that, for all sufficiently large n, with probability at least 1 − 2n−4 , simultaneously for all s ∈ {0, 1, 2}, 

sup

b ∇s ℓ(x) − ∇s ℓ(x)

≤ Cj hβ−s +

x∈B(µj ,rj )



log n  . nhd+2s

Lemma 9 (High-probability analytic local event). Under Assumptions 2, 3, and 1, for each j ∈ [M ], Pr Alocal,j ≥ 1 − 4n−4 n for all sufficiently large n. The next steps use Alocal,j to control the stopped local dynamics and then intersect it with n stay,j local,j En,T to recover the full event Gn,T . Lemma 10 (Local score-contribution bound). Let us write pmin := minj∈[M ] pmin,j , and C∗ := 2GK −(d+1) . For pmin h gi (x) 1 x − Xi qi (x) := , gi (x) := d+1 ∇K , pb(x) h h one has, for every j ∈ [M ], on Alocal,j , n sup x∈B(µj ,rj )

max ∥qi (x)∥ ≤ C∗

1≤i≤n

for all sufficiently large n. In particular, C∗ ≍ h−(d+1) , which is the local scale used to control the privacy perturbation in the proof of Theorem 4.3. 36

Proposition A.1 (Local inactivity of stabilization). Suppose β > 2 and let Assumptions 1, 2, and 3 hold. Let us fix pmin := minj∈[M ] pmin,j ,, 0 < cfloor ≤ 12 , and set pfloor = cfloor pmin . Let us choose A > 0 so that A > 2 maxj∈[M ] supx∈B(µj ,rj ) ∥∇p(x)∥. Then, for each j ∈ [M ], on Alocal,j and for n all sufficiently large n, pb(x) ≥ pfloor ,

∥∇pb(x)∥ ≤ A,

Consequently, sbA,pfloor ;X (x) = ∇ log pb(x)

x ∈ B(µj , rj ).

for all x ∈ B(µj , rj ).

Proof of Proposition A.1. Fix j ∈ [M ] and work on Alocal,j . By Definition 5, n pb(x) ≥ cj =

pmin,j pmin ≥ ≥ pfloor 2 2

for all x ∈ B(µj , rj ),

so the denominator floor is inactive on the whole basin. It remains to show that gradient clipping is also inactive. Let an :=

b |ℓ(x) − ℓ(x)|,

sup

bn :=

x∈B(µj ,rj )

b ∥∇ℓ(x) − ∇ℓ(x)∥.

sup x∈B(µj ,rj )

On Alocal,j , we have an → 0 and bn → 0. Since pb(x) = p(x)ebℓ(x)−ℓ(x) , n |pb(x) − p(x)| ≤ pmax,j ean − 1 → 0,

sup x∈B(µj ,rj )

where pmax,j := supx∈B(µj ,rj ) p(x) < ∞. Also, using b ∇pb(x) = pb(x)∇ℓ(x),

∇p(x) = p(x)∇ℓ(x),

we obtain b b ∥∇pb(x) − ∇p(x)∥ ≤ |pb(x) − p(x)| ∥∇ℓ(x)∥ + p(x) ∥∇ℓ(x) − ∇ℓ(x)∥

≤ |pb(x) − p(x)| ∥∇ℓ(x)∥ + bn + pmax,j bn . Taking suprema over x ∈ B(µj , rj ), and using compactness of B(µj , rj ) together with continuity of p and ℓ, gives sup ∥∇pb(x) − ∇p(x)∥ → 0 on Alocal,j . n x∈B(µj ,rj )

This convergence is the s = 1 derivative control in Lemma 6, used under Assumption 3. Therefore sup x∈B(µj ,rj )

∥∇pb(x)∥ ≤

∥∇p(x)∥ + o(1) < A

sup x∈B(µj ,rj )

for all sufficiently large n, by the choice of A. Hence, on Alocal,j and for all sufficiently large n, clipA (∇pb(x)) = ∇pb(x) and max{pb(x), pfloor } = n pb(x) for all x ∈ B(µj , rj ). Therefore sbA,pfloor ;X (x) = ∇ log pb(x) as claimed. 37

for all x ∈ B(µj , rj ),

We now collect the local bounds used to control the stabilized minibatch field. Lemma 10 gives a uniform bound on the samplewise score contributions, Lemma 11 gives a conditional secondmoment bound for the stabilized minibatch field, and Lemma 12 upgrades this to a high-probability fluctuation bound along adapted trajectories. Lemma 11 (Second-moment control of the minibatch stabilized field). Assume β > 2 and suppose Assumptions 1, 2, and 3 hold. Fix j ∈ [M ] and work on Alocal,j . Then there exists a constant n Cj,ζ > 0 such that, for every deterministic x ∈ B(µj , rj ), i

≤ E ∥sbA,pfloor ;Bt (x) − ∇ log pb(x)∥2 x, X , Alocal,j n

Cj,ζ . mhd+2

Lemma 12 (High-probability minibatch fluctuation of the stabilized field). Assume β > 2 and suppose Assumptions 1, 2, and 3 hold. Fix j ∈ [M ], assume mhd+2 ≥ C, log(eT n) −1 be any process adapted to the algorithmic history for a sufficiently large constant C > 0. Let (xt )Tt=0 such that xt ∈ B(µj , rj ) for all t. Then there exists a constant Cj,mb > 0 such that



Pr max ∥sbA,pfloor ;Bt (xt ) − ∇ log pb(xt )∥ ≤ Cj,mb 0≤t≤T −1



log(eT n)  ≥ 1 − T (eT n)−6 X , Alocal,j n mhd+2

for all sufficiently large n. Lemma 13 (Local score bounds for kn ). Assume β > 2 and suppose Assumptions 1, 2, 3, 4, and 4 hold. Fix j ∈ [M ]. For xt ∈ B(µj , rj ), write δt := xt − µj ,

b kn (x) := ∇ log pb(x) = ∇ℓ(x).

, there exists a constant Cj > 0 such that for all sufficiently large n, Then, on Alocal,j n 3αj log n ⟨δt , kn (xt )⟩ ≤ − ∥δt ∥2 + Cj h2(β−1) + + C hβ−2 + j 8 nhd+2

and

∥kn (xt )∥2 ≤ 2Hj2 ∥δt ∥2 + Cj h2(β−2) +

log n ∥δt ∥2 , nhd+4

log n log n 2 2(β−1) ∥δ ∥ + C h + . t j nhd+4 nhd+2

Lemma 14 (Finite-population second moment for the local score average). Assume the hypotheses of Lemma 13. For 1 X gi (xt ) q̄t := qi (xt ), qi (xt ) = , m i∈B pb(xt ) t

there exist constants C3,j , C4,j > 0 such that, on Alocal,j , for every xt ∈ B(µj , rj ) and all sufficiently n large n, log n 1 2 2(β−1) E[∥q̄t ∥2 | xt , X , Alocal,j ] ≤ C ∥δ ∥ + C h + + . 3,j t 4,j n nhd+2 mhd+2 38

Lemma 15 (Minibatch remainder bound). Let us define rt := sbA,pfloor ;Bt (xt ) − q̄t ,

q̄t :=

1 X qi (xt ). m i∈B t

Then, on Alocal,j and for xt ∈ B(µj , rj ), n E[∥rt ∥2 | xt , X , Alocal,j ]≲ n

1 . mhd+2

For the local proofs, fix j ∈ [M ] and write ℓ = log p. Define Hj :=

∥∇2 ℓ(x)∥ < ∞,

sup

Lj :=

x∈B(µj ,rj )

∥∇ log p(x)∥.

sup x∈B(µj ,rj )

The finiteness of Hj follows from Lemma 4 and compactness of B(µj , rj ). Fix any deterministic Bj > Lj , and set η0,j :=

αj rj2 , 8Bj2

η1,j :=

αj , 48(Hj2 + 1)

η̄j := min{η0,j , η1,j }.

The next proposition gives the fixed-bandwidth local error bound on the stay-in-basin event, conditional on the static analytic event Alocal,j . Basin retention is handled afterward and will then n allow us to pass from the stopped estimate to the full local theorem. Proposition A.2 (General stopped local error bound). Assume β > 2 and suppose Assumptions 1, 2, 3, 4, and 4 hold. Let x0 ∈ B(µj , rj ) for some j ∈ [M ], choose pfloor and A as in Proposition A.1, and let 0 < η ≤ η1,j . Then there exist constants κj > 0 and Cj > 0, independent of n, such that for all sufficiently large n,

E ∥xT − µj ∥ 1E stay,j n,T

X , Alocal,j n

2(β−1)

≤ (1−κj η) ∥x0 −µj ∥ +Cj h

log n 1 + + + ηdσ 2 . nhd+2 mhd+2

throughProof of Proposition A.2. Condition on the observed sample X and on the event Alocal,j n out. For t = 0, . . . , T , let us define stay,j Et := En,t = {xs ∈ B(µj , rj ) for all s = 0, . . . , t}.

and the filtration Ft := σ(B0 , Z0 , . . . , Bt−1 , Zt−1 , x0 ),

t ≥ 0.

Then Et ∈ Ft . The proof tracks the squared error only up to the first exit from the local basin. The indicators 1Et therefore localize the recursion to the regime in which the basinwise analytic controls from the preceding lemmas are valid. We write δt := xt − µj ,

q̄t :=

1 X qi (xt ), m i∈B t

rt := sbA,pfloor ;Bt (xt ) − q̄t ,

so that zt ∼ N (0, σ 2 Id ).

xt+1 = xt + η(q̄t + rt + zt ), Since Et+1 ⊆ Et ,

∥δt+1 ∥2 1Et+1 ≤ 1Et δt + η(q̄t + rt + zt ) . Taking conditional expectation given Ft , and writing D

T2 := δt , E[q̄t | Ft , X , Alocal,j ] , n therefore yields

]⟩ ≤ 1Et ∥δt ∥2 + 2ηT2 + 2η⟨δt , E[rt | Ft , X , Alocal,j E ∥δt+1 ∥2 1Et+1 Ft , X , Alocal,j n n

+ η 2 E[∥q̄t + rt + zt ∥2 | Ft , X , Alocal,j ] . n We next separate the one-step recursion into its principal drift term, the stabilization remainder, and the quadratic second-moment term. On Et , one has xt ∈ B(µj , rj ), so the local lemmas apply. First, E[q̄t | Ft , X , Alocal,j ] = ∇ log pb(xt ) =: kn (xt ), n and Lemma 13 gives log n 3αj ∥δt ∥2 + Cj h2(β−1) + 8 nhd+2

T2 = ⟨δt , kn (xt )⟩ ≤ −



+ Cj hβ−2 +



log n  ∥δt ∥2 . nhd+4

Next, i

, E[rt | Ft , X , Alocal,j ] = E sbA,pfloor ;Bt (xt ) − ∇ log pb(xt ) Ft , X , Alocal,j n n so Lemma 11 and Jensen’s inequality yield E[rt | Ft , X , Alocal,j ] n

≤ E[∥rt ∥ | Ft , X , Alocal,j ]≤ n

E[∥rt ∥

| Ft , X , Alocal,j ] ≤ Cj n

1 . mhd+2

Hence, by Young’s inequality, D

2η δt , E[rt | Ft , X , Alocal,j ] ≤ n

αj η 1 ∥δt ∥2 + Cj η . 8 mhd+2

It remains to control the quadratic contribution. Here the average score term and the stabilization remainder are handled by the second-moment bounds from Lemmas 14 and 15, while the Gaussian perturbation contributes the explicit dσ 2 term. For the quadratic term, by (a + b + c)2 ≤ 3a2 + 3b2 + 3c2 , E[∥q̄t + rt + zt ∥2 | Ft , X , Alocal,j ] ≤ 3 E[∥q̄t ∥2 | Ft , X , Alocal,j ] + 3 E[∥rt ∥2 | Ft , X , Alocal,j ] + 3dσ 2 . n n n

By Lemmas 14 and 15, there exist constants C3,j , C4,j > 0 such that, on Et , E[∥q̄t ∥2 | Ft , X , Alocal,j ] + E[∥rt ∥2 | Ft , X , Alocal,j ] n n

≤ C3,j ∥δt ∥2 + C4,j h2(β−1) +

log n 1 + d+2 nh mhd+2

for all sufficiently large n. Combining the previous displays and using Assumption 3, which gives s

hβ−2 +

log n = o(1), nhd+4

together with 0 < η ≤ η1,j , we obtain constants κj > 0 and Cj′ > 0 such that h

E ∥δt+1 ∥2 1Et+1 Ft , X , Alocal,j ≤ (1 − κj η) ∥δt ∥2 1Et n + Cj′ η

2(β−1)

1 log n + + nhd+2 mhd+2

+ Cj′ η 2 dσ 2

for all sufficiently large n. Let us now define h

2 local,j ∆tr . t := E ∥δt ∥ 1Et X , An

Taking expectations in the last display gives tr ′ ∆tr t+1 ≤ (1 − κj η)∆t + Cj η

2(β−1)

log n 1 + + d+2 nh mhd+2

+ Cj′ η 2 dσ 2 .

Unrolling the recursion yields T 2 ′′ ∆tr T ≤ (1 − κj η) ∥x0 − µj ∥ + Cj

where Cj′′

2(β−1)

0 is independent of n, m, h, T . Since ∆tr T

log n 1 + + + ηdσ 2 , d+2 nh mhd+2

E ∥xT − µj ∥2 1E stay,j n,T

X , Alocal,j n

, this

proves the proposition. The next two results control the perturbation term uniformly up to exit and then show that, under the stated tuning conditions, the trajectory remains in B(µj , rj ) with high probability. Lemma 16 (Uniform perturbation control up to exit). Assume β > 2. Suppose Assumptions 1, 2, 3, 4, and 4 hold, and assume mhd+2 → ∞. log(eT n) −1 Fix j ∈ [M ], and let (xt )Tt=0 be any process adapted to the algorithmic filtration such that x0 ∈ B(µj , rj ). For each t = 0, . . . , T − 1, let us define (j)

:= {xs ∈ B(µj , rj ) for all s = 0, . . . , t},

(j)

et x

:= xt 1E (j) + x0 1(E (j) )c . t

Also recall C∗ := (2GK /pmin ) h−(d+1) from Lemma 10. Define 

(j) Ξn,T,m,h := C1,j hβ−1 +

log n + nhd+2



C∗ T d polylog(T, n, δ) log(eT n)  + C2,j , d+2 mh nεmodes

where C1,j , C2,j > 0 are deterministic constants independent of n, m, h, T . Then, on Alocal,j , n

max

0≤t≤T −1

(j) (j) e t ) + zt et ) − ∇ log p(x sbA,pfloor ;Bt (x

(j) ≤ Ξn,T,m,h

X , Alocal,j n

≥ 1 − 2T (eT n)−6 .

Proposition A.3 (High-probability basin retention). Assume β > 2. Suppose Assumptions 1, 2, mhd+2 3, 4, and 4 hold, and assume log(eT n) → ∞. Let x0 ∈ B(µj , rj ) for some j ∈ [M ], and choose pfloor and A as in Proposition A.1. Define stay,j := {xt ∈ B(µj , rj ) for all t = 0, . . . , T }. En,T (j)

With Bj and η0,j as defined in the local proof notation above, if 0 < η ≤ η0,j and Ξn,T,m,h ≤ α r , for all sufficiently large n min Bj , j4 j , then on Alocal,j n

stay,j X , Alocal,j ≥ 1 − 2T (eT n)−6 . Pr En,T n

Proof of Proposition A.3. We condition on X and on Alocal,j throughout. For t = 0, . . . , T , define n (j)

:= {xs ∈ B(µj , rj ) for all s = 0, . . . , t}.

For t = 0, . . . , T − 1, define the stopped surrogate (j)

et x

Let (j) Fn,T :=

max

0≤t≤T −1

:= xt 1E (j) + x0 1(E (j) )c . t

(j) (j) et ) − ∇ log p(x e t ) + zt sbA,pfloor ;Bt (x

(j) ≤ Ξn,T,m,h

By Lemma 16,

(j)

Pr Fn,T X , Alocal,j ≥ 1 − 2T (eT n)−6 . n (j)

We show that on Fn,T , the actual trajectory never leaves B(µj , rj ). We proceed by induction (j)

(j)

on t. The base case E0 holds because x0 ∈ B(µj , rj ) by assumption. Now suppose Et some t ∈ {0, . . . , T − 1}. Then (j) et = xt . x (j)

(j)

On Fn,T , there therefore exists a vector ut with ∥ut ∥ ≤ Ξn,T,m,h such that sbA,pfloor ;Bt (xt ) + zt = ∇ log p(xt ) + ut .

holds for

By the algorithm update,

xt+1 = xt + η sbA,pfloor ;Bt (xt ) + zt = xt + η ∇ log p(xt ) + ut . We write δt := xt − µj . Then ∥xt+1 − µj ∥2 = ∥δt ∥2 + 2η⟨δt , ∇ log p(xt )⟩ + 2η⟨δt , ut ⟩ + η 2 ∥∇ log p(xt ) + ut ∥2 . By Lemma 5, ⟨δt , ∇ log p(xt )⟩ ≤ −

αj ∥δt ∥2 . 2

Also, (j)

⟨δt , ut ⟩ ≤ ∥δt ∥ ∥ut ∥ ≤ rj Ξn,T,m,h , (j)

and, since ∥∇ log p(xt )∥ ≤ Lj < Bj and Ξn,T,m,h ≤ Bj , ∥∇ log p(xt ) + ut ∥2 ≤ (∥∇ log p(xt )∥ + ∥ut ∥)2 ≤ 4Bj2 . Therefore, (j)

∥xt+1 − µj ∥2 ≤ (1 − αj η)∥δt ∥2 + 2ηrj Ξn,T,m,h + 4η 2 Bj2 . (j)

Because Et

holds, ∥δt ∥ ≤ rj . Using αj rj 4

and

2ηrj Ξn,T,m,h ≤

αj ηrj2 , 2

(j)

Ξn,T,m,h ≤ we obtain (j)

Hence ∥xt+1 − µj ∥2 ≤ (1 − αj η)rj2 +

η ≤ η0,j =

4η 2 Bj2 ≤

αj rj2 , 8Bj2

αj ηrj2 . 2

αj ηrj2 αj ηrj2 + = rj2 . 2 2

(j)

Thus xt+1 ∈ B(µj , rj ), so Et+1 holds. (j) (j) (j) stay,j By induction, Et holds for all t = 0, . . . , T on Fn,T . Therefore Fn,T ⊆ En,T , and so

stay,j Pr En,T X , Alocal,j ≥ 1 − 2T (eT n)−6 . n

This proves the proposition. We now combine the stopped local recursion with basin retention to derive the stated rate and identify the optimizing bandwidth.

local,j stay,j Proof of Theorem 4.3. Recall from Definition 5 that Gn,T = Alocal,j ∩En,T . By the notation used n in the proof of the local convergence, η̄j = min{η0,j , η1,j }. Hence 0 < η ≤ η̄j implies 0 < η ≤ η0,j and 0 < η ≤ η1,j . By Lemma 9, Pr Alocal,j ≥ 1 − 4n−4 n

for all sufficiently large n. We first derive the local error bound for a generic bandwidth h, and only then optimize in h. Conditional on (X , Alocal,j ), Proposition A.1 shows that the stabilized score agrees with the ordinary n log-KDE score throughout B(µj , rj ) for all sufficiently large n. On the same event, Lemma 10 gives the uniform samplewise bound max ∥qi (x)∥ ≤ C∗ ,

sup x∈B(µj ,rj )

C∗ =

1≤i≤n

2GK −(d+1) h . pmin

This is the local quantity that will control the privacy contribution. By Proposition A.2, there exists a constant Cj > 0, independent of n, such that on Alocal,j , n

E ∥xT − µj ∥ 1E stay,j n,T

X , Alocal,j n

2(β−1)

≤ (1−κj η) ∥x0 −µj ∥ +Cj h

log n 1 + + + ηdσ 2 . d+2 nh mhd+2

We next bound the privacy term. Since pfloor and A are fixed as in Proposition A.1, the definition (8) gives 1/2

√ ∆h,corr (A, pfloor ) = 2 2

1/2

I1 AI h−(d+1) + 2 0 h−d . pfloor pfloor

Because h → 0 under Assumption 3, one has h ≤ 1 for all sufficiently large n, hence h−d ≤ h−(d+1) . Therefore ∆h,corr (A, pfloor ) ≲ h−(d+1) ≍ C∗ . Using (9), the definition of εiter , and the fact that η is fixed independently of n, we obtain ηdσ 2 ≲

T d polylog(n, δ) −2(d+1) h . n2 ε2modes

Also, since m ≍ n/ log n, 1 log n ≲ . d+2 mh nhd+2 Thus, after absorbing constants,

E ∥xT − µj ∥2 1E stay,j X , Alocal,j ≤ (1 − κj η)T ∥x0 − µj ∥2 + Cj Rn (h), n

(16)

n,T

where Rn (h) := h2(β−1) +

log n + An h−2(d+1) , nhd+2

An :=

T d polylog(n, δ) . n2 ε2modes

We now optimize Rn (h). We first balance the approximation term and the nonprivate stochastic term: log n h2(β−1) ≍ . nhd+2 44

This gives hnp ≍

log n

1 d+2β

2(β−1)

and

Thus Rn (hnp ) ≲

h2(β−1) ≍ np

log n 2(β−1) d+2β

n The nonprivate choice remains optimal whenever An h−2(d+1) ≲ np

log n log n d+2β . ≍ d+2 n nhnp

+ An h−2(d+1) . np

log n 2(β−1) d+2β

which is equivalent, up to logarithmic factors and using T = CT log n, to εmodes ≳ εthr . In that regime, log n 2(β−1) d+2β Rn (hnp ) ≲ . n Let us next balance the approximation term and the privacy term: h2(β−1) ≍ An h−2(d+1) . This gives 1 T d polylog(n, δ) 1 2d+2β , hdp ≍ An2d+2β = 2 2 n εmodes and hence 2(β−1)

hdp

−2(d+1)

≍ An hdp

β−1

≍ And+β =

T d polylog(n, δ) β−1 d+β

n2 ε2modes

If εmodes ≲ εthr , then the remaining middle term is dominated by the common order of the first and third terms, so T d polylog(n, δ) β−1 d+β Rn (hdp ) ≲ . 2 2 n εmodes Combining the two regimes, for the bandwidth choice h = hopt in (11), Rn (hopt ) ≲

log n 2(β−1) d+2β

T d polylog(n, δ) β−1 d+β

n2 ε2modes

(17)

We next verify the stay-in-basin hypotheses of Proposition A.3. Since m ≍ n/ log n and T = CT log n, the first requirement mhd+2 →∞ log(eT n)

For h = hnp ≍ log n/n

1 d+2β

is equivalent to

nhd+2 → ∞. (log n)2

, this gives 2(β−1) d+2 nhd+2 np −2+ d+2β → ∞. ≍ n d+2β (log n) 2 (log n)

For

h = hdp ≍ An2d+2β ,

An :=

T d polylog(n, δ) , n2 ε2modes

the privacy-dominated condition εmodes ≲ εthr is equivalent, up to logarithmic factors, to An ≳

log n 2(d+β) d+2β

Therefore d+2 2d+2β

nhd+2 dp

nAn = ≳ (log n)2 (log n)2

n logn n

d+2

d+2β 2(β−1)

= n d+2β (log n)

(log n)2

d+2 −2+ d+2β

→ ∞.

Thus the minibatch-growth condition from Proposition A.3 holds in both regimes. It remains to verify that αj rj (j) Ξn,T,m,h ≤ min Bj , 4 for all sufficiently large n. Since 

(j) Ξn,T,m,h = C1,j hβ−1 +

log n + nhd+2



C∗ T d polylog(T, n, δ) log(eT n)  + C2,j , mhd+2 nεmodes

it is enough to show that each term vanishes in both regimes. For h = hnp , hβ−1 np =

log n β−1

d+2β

→ 0,

β−1

log n log n d+2β → 0, = n nhd+2 np

and, since m ≍ n/ log n and log(eT n) ≍ log n, s

−(d+1)

Also, C∗ ≍ hnp

v u

2 d+2 β−1 log(eT n) u t log n = n− d+2β (log n)1− 2(d+2β) → 0. ≲ d+2 d+2 mhnp nhnp

, so p

1/2 C∗ T d polylog(T, n, δ) p ≍ An h−(d+1) = An h−2(d+1) . np np nεmodes

In the nonprivate regime, An h−2(d+1) ≲ np

log n 2(β−1) d+2β

β−1

hence

C∗ T d polylog(T, n, δ) log n d+2β → 0. ≲ nεmodes n

β−1

2d+2β For h = hdp , hβ−1 and, since hdp is defined by balancing approximation and privacy, dp = An

−(d+1)

An hdp

−2(d+1) 1/2

= An hdp

≍ hβ−1 dp .

Moreover, the minibatch-growth condition already proved implies v u u log(eT n) t → 0,

mhd+2 dp

and therefore also

v u u log n t → 0. d+2

nhdp

(j)

Thus every term in Ξn,T,m,h vanishes in both regimes. Therefore Ξn,T,m,h → 0, and so for all sufficiently large n, αj rj (j) . Ξn,T,m,h ≤ min Bj , 4 Applying Proposition A.3 with the corresponding bandwidth choice now yields

stay,j Pr En,T X , Alocal,j ≥ 1 − n−5 n

for all sufficiently large n. Therefore local,j stay,j Pr Gn,T = Pr Alocal,j ∩ En,T ≥ 1 − 4n−4 − n−5 ≥ 1 − 5n−4 . n

local,j stay,j Since Gn,T = Alocal,j ∩ En,T , we have n

local,j ≤ E ∥xT − µj ∥2 X , Gn,T

E ∥xT − µj ∥2 1E stay,j X , Alocal,j n

n,T

stay,j En,T

X , Alocal,j n

Using the lower bound 1 − n−5 in the denominator and the bound (16) together with (17), we obtain h i (1 − κj η)T ∥x0 − µj ∥2 + Cj Rn (hopt ) local,j E ∥xT − µj ∥2 X , Gn,T ≤ . 1 − n−5 Finally, since T = CT log n, (1 − κj η)T ≤ e−κj ηT = n−κj ηCT . ⋆ > 0 so that We choose CT,j

2(β − 1) ⋆ κj ηCT,j > d + 2β

(1 − κj η) = o

and hence

log n 2(β−1)

d+2β

⋆ . Absorbing this transient term and the factor (1 − n−5 )−1 into the constants for every CT ≥ CT,j completes the proof:

E ∥xT − µj ∥

A.3

local,j X , Gn,T

≤ CnonDP,j

log n 2(β−1) d+2β

+ CDP,j

T d polylog(n, δ) β−1 d+β

n2 ε2modes

Global convergence: DAP coverage, basinwise trajectories, and post-processing

We next turn from basinwise control to global recovery. The first group of lemmas analyzes the DAP design and shows that the initialization pool covers all modal basins with high probability. The final step combines this coverage event with the local theorem and the deterministic post-processing conditions used by the merge rule to prove Theorem 4.5. Throughout the DAP coverage argument, for j ∈ [M ] and zr ∈ Zn , write Gj,n := Zn ∩ B(µj , ρinit /4),

ur :=

n 1X 1{∥Xi − zr ∥ ≤ hDAP }, n i=1

mr := Pr(∥X − zr ∥ ≤ hDAP ).

Lemma 17 (Geometric consequences of the DAP design). Assume Assumption 2, and work under the DAP design specified above. Then, for all sufficiently large n, the following hold: (i) Gj,n ̸= ∅ for every j ∈ [M ]. (ii) If a ∈ B(µj , ρinit /4), then B(µj , ρinit /4) ⊆ B(a, ρinit ). (iii) If a ∈ B(µj , ρinit /4), then B(a, ρinit ) ∩ B(µi , ρinit /4) = ∅

for every i ̸= j.

Lemma 18 (Within-basin local-mass ordering). Suppose Assumptions 2, 4 hold. Work under the DAP design specified above. Then, for each j ∈ [M ], there exist constants ∆j > 0 and nj ∈ N such that, for all n ≥ nj , inf r: zr ∈Zn ∩B(µj ,ρinit /4)

mr ≥

ms + ∆j hdDAP ρ2init .

sup s: zs ∈Zn ∩(B(µj ,rj )\B(µj ,ρinit /4))

Lemma 19 (Uniform concentration of the DAP utilities). Suppose Assumption 2 holds, and work under the DAP design specified above. Fix cconc > 0. Then there exists CDAP,tail > 0, independent of εinit , such that, for all sufficiently large n,

max

1≤r≤Ncand

|ur − mr | ≤ cconc hdDAP ρ2init

≥ 1 − CDAP,tail n−4 .

Lemma 20 (DAP competitive-region complexity). Suppose Assumptions 2, 4 hold. Work under the DAP design specified above. Fix m ∈ [M ]. For a sufficiently small fixed γ0 > 0, define

Cm,n := zr ∈ Zn \ B(µm , rm ) : mr ≥

inf

s: zs ∈Gm,n

ms − γ0 hdDAP ρ2init .

Then, for all sufficiently large n, [

Cm,n ⊆

B(µj , rj ),

j∈[M ]\{m}: p(µj )≥p(µm )

and Cm,n can be covered by at most Lcap M log n Euclidean balls of radius ρinit /2, where Lcap > 0 is uniform over m ∈ [M ] and independent of n and εinit . Proof of Proposition 4.4. Let cconc ≤ γ0 /8, where γ0 is the constant in Lemma 20. Define

Econc :=

max

1≤r≤Ncand

|ur − mr | ≤ cconc hdDAP ρ2init

Lemma 19 gives c Pr(Econc ) ≤ CDAP,tail n−4 .

Work on Econc , and put

nεinit . 2k For any active set Aℓ , any nonempty T ⊆ Aℓ , and any S ⊆ Aℓ , the exponential mechanism gives λ :=

exp(λur ) |S| exp{λ supr: zr ∈S ur } P ≤ q: zq ∈Aℓ exp(λuq ) q: zq ∈T exp(λuq )

Pr(Jℓ ∈ S | Aℓ ) = P r: zr ∈S "

(

≤ |S| exp −λ

inf uq − sup ur

q: zq ∈T

r: zr ∈S

For each j ∈ [M ], the separation argument in Lemma 20 gives a constant ηj > 0 such that p(x) ≤ p(µj ) − ηj .

sup x∈Q\∪M i=1 B(µi ,ri )

Choose sj ∈ (0, rj /3) so that inf

∥x−µj ∥≤sj

p(x) ≥ p(µj ) − ηj /4.

The local-mass expansion gives a constant Clm > 0 such that, for all zr ∈ Zn ∩ Q, mr − hdDAP Vol(B(0, 1))p(zr ) ≤ Clm hd+1 DAP . By the uniform DAP separation constants, η⋆ := inf 1≤j≤M ηj > 0. Set cout :=

1 Vol(B(0, 1))η⋆ . 2

Whenever 2Clm hDAP ≤

1 Vol(B(0, 1))η⋆ , 4

we have inf r: zr ∈Zn ∩∪M j=1 B(µj ,sj )

mr ≥

sup

ms + cout hdDAP .

s: zs ∈Zn \∪M j=1 B(µj ,rj )

For each j, take a 4ρinit -separated subset of B(µj , sj /2) of cardinality at least cρ−d init , and project each point to its nearest grid point in Zn . Choose n0 so that, for all n ≥ n0 , √ d ρinit ρinit 1 c0 hDAP ≤ , < min sj , ρinit < . 2 8 8 2 1≤j≤M 3 The projected points lie in B(µj , sj ). Points projected from the same mode remain 2ρinit -separated because ρinit > 2ρinit . ∥z − z ′ ∥ ≥ 4ρinit − 2 · 8 If z and z ′ are projected from different modes i ̸= j, then 2 2 ∥z − z ′ ∥ ≥ ∥µi − µj ∥ − si − sj ≥ ∥µi − µj ∥ ≥ c0 > 2ρinit . 3 3

Let Pn be the union of the projected grid points over j ∈ [M ]. Since ρinit ≍ (log n)−1/d , there is cpack > 0 such that |Pn | ≥ cpack M ρ−d init ≥ cpack M log n. No suppression ball of radius ρinit can remove two points of Pn . The upper implicit constant in k ≍ M log n is chosen so that cpack M log n. k≤ 2 Thus, before every draw ℓ ≤ k, at least one point of Pn remains active. Set Γ := min{cout /2, 3γ0 /4}. We next use the lower bound on εinit explicitly. The notation εinit ≳ M n−2β/(d+2β) polylog(n) is used here in the following sufficient form: there exists a constant Cε,init > 0 such that εinit ≥ Cε,init M n−2β/(d+2β) (log n)2−d/(d+2β)+2/d . Since k ≍ M log n, there is Ck > 0 such that k ≤ Ck M log n. Since hDAP ≍ (log n/n)1/(d+2β) and ρinit ≍ (log n)−1/d , there are ch , cρ > 0 such that hdDAP ≥ ch Therefore

log n n

d/(d+2β)

ρ2init ≥ cρ (log n)−2/d .

nεinit d h ρ2 2k DAP init n ≥ · Cε,init M n−2β/(d+2β) (log n)2−d/(d+2β)+2/d 2Ck M log n

λhdDAP ρ2init =

× ch =

log n n

d/(d+2β)

cρ (log n)−2/d

Cε,init ch cρ log n. 2Ck

Choose Cε,init ≥

18Ck . Γch cρ

Then ΓλhdDAP ρ2init ≥ 9 log n. Fix any draw ℓ ≤ k. Let T be any active singleton contained in Pn , and set S = Aℓ \

M [ j=1

B(µj , rj ).

On Econc , the modal-comparator gap gives inf uq − sup ur ≥ inf mq − sup mr − 2cconc hdDAP ρ2init

q: zq ∈T

r: zr ∈S

≥ cout hdDAP − 2cconc hdDAP ρ2init . Increase n0 , if necessary, so that 2cconc ρ2init ≤ cout /2 and ρ2init ≤ 1 for all n ≥ n0 . Then inf uq − sup ur ≥ ΓhdDAP ρ2init .

q: zq ∈T

r: zr ∈S

Since |S| ≤ Ncand ≤ Ch−d DAP ≤ Cn, the exponential-mechanism bound gives 

PrJℓ ∈ /

M [



B(µj , rj ) | Aℓ  ≤ Cn exp −ΓλhdDAP ρ2init ≤ Cn−8 .

j=1

A union bound over k ≍ M log n draws gives 

/ Pr∃ℓ ≤ k : aℓ ∈

M [



B(µj , rj ) Econc  ≤ CM (log n)n−8 ≤ Cn−7 ,

j=1

where the last inequality uses the implicit growth regime for M in the DAP design. Let Eloc denote S the event that all selected anchors lie in M j=1 B(µj , rj ). It remains to show that every modal basin is reached. Let λ(1) > λ(2) > · · · > λ(R) be the distinct values among {p(µj ) : j ∈ [M ]}, and define Hr := {j ∈ [M ] : p(µj ) = λ(r) },

r ∈ [R].

We work on Econc ∩ Eloc . By Lemma 20, M [

Cm,n ⊆

m=1

M [

B(µj , rj ) ⊆ Q.

j=1

Since the public box Q is fixed and ρinit ≍ (log n)−1/d , this union of competitive sets can be covered by at most Lpeel M log n Euclidean balls of radius ρinit /2, with Lpeel independent of M , n, and εinit . The constants implicit in the DAP design choice k ≍ M log n are fixed so that, for all n ≥ n0 , Lpeel M log n + M ≤ k ≤

cpack M log n. 2

The upper bound is the packing requirement used above to keep an active modal comparator throughout the k draws; the lower bound is the peeling requirement used below. We prove coverage level by level. Suppose that all modes in H1 , . . . , Hr−1 have been hit. Let Vr ⊆ Hr be the modes in the current level that have already been hit, and set Ur := Hr \ Vr . If Ur = ∅, move to level r + 1. Otherwise define TUr ,n :=

[

Gj,n ,

BUr :=

[ j∈Ur

j∈Ur

B(µj , rj ).

For every j ∈ Ur , Gj,n remains active until mode j is hit. Indeed, on Eloc , every earlier anchor lies in some B(µi , ri ). If i ̸= j, a ∈ B(µi , ri ), and z ∈ Gj,n , then ∥z − a∥ ≥ ∥µi − µj ∥ − ∥a − µi ∥ − ∥z − µj ∥ ≥

c0 ρinit − . 2 4

Increasing n0 , if necessary, so that ρinit ≤ c0 /5, the last display is larger than ρinit . Thus suppressions from other basins do not remove Gj,n , and no anchor has yet fallen in B(µj , rj ) because j ∈ Ur . Choose m⋆r ∈ Ur such that inf

q: zq ∈Gm⋆ ,n

mq =

inf

q: zq ∈TUr ,n

mq .

For a draw made before level r is completed, set S0 := Aℓ \ (BUr ∪ Cm⋆r ,n ). Then the definition of Cm⋆r ,n gives sup ms < s: zs ∈S0

inf

q: zq ∈TUr ,n

On Econc , inf

q: zq ∈TUr ,n

uq − sup us ≥ s: zs ∈S0

mq − γ0 hdDAP ρ2init .

3γ0 d h ρ2 ≥ ΓhdDAP ρ2init . 4 DAP init

The exponential-mechanism bound therefore gives

Pr(Jℓ ∈ S0 | Aℓ ) ≤ Cn exp −ΓλhdDAP ρ2init ≤ Cn−8 . Thus, except on an event of conditional probability at most Cn−8 , each draw before the current level is completed either lands in BUr , which hits at least one mode in Ur , or lands in the active part of Cm⋆r ,n . By Lemma 20, the latter set lies in basin neighborhoods whose heights are at least λ(r) . Since points in Aℓ \ BUr are outside the unvisited level-r basins, any such active competitive point lies in a higher-level basin or in an already-hit same-level basin. If a selected anchor aℓ lies in one of the fixed competitive-cover balls with center y, then every grid point z in that same ball is removed at the next suppression step, since ∥z − aℓ ∥ ≤ ∥z − y∥ + ∥aℓ − y∥ ≤

ρinit ρinit + = ρinit . 2 2

Because the packing argument above keeps the active set nonempty throughout the k draws, the reset clause in Algorithm 2 is not invoked on this event. Therefore a competitive-cover ball, once suppressed, is not charged again. Over the full run, at most Lpeel M log n non-hit draws can be charged to active competitive-cover balls. Since at most M successful hits are needed to hit all basins and k ≥ Lpeel M log n + M , failure to hit every basin by time k implies that at least one exceptional event Jℓ ∈ S0 occurred. A union bound over k ≍ M log n draws yields

Pr ∃j ∈ [M ] : I ∩ B(µj , rj ) = ∅ Econc ∩ Eloc ≤ CM (log n)n−8 ≤ Cn−7 , 52

again using the implicit growth regime for M in the DAP design. Combining the concentration, localization, and coverage bounds gives 

PrI ⊆



M [

B(µj , rj ),

I ∩ B(µj , rj ) ̸= ∅ for every j ∈ [M ]

j=1

≥ 1 − CDAP,tail n−4 − Cn−7 − Cn−7 ≥ 1 − Cinit,cov n−2 .

Lemma 21 (High-probability static global event). Under Assumptions 1, 2, and 3, Pr Aglobal ≥ 1 − Cglobal,stat n−2 n

for all sufficiently large n, for some constant Cglobal,stat > 0. Lemma 22 (Post-processing preserves basinwise rates). Recall Ij := {ℓ ∈ [k] : x0,ℓ ∈ B(µj , rj )} for j ∈ [M ]. Fix an event E and deterministic bounds Rn,j such that h

E ∥xT,ℓ − µj ∥2 X , E ≤ Rn,j

for every j ∈ [M ] and every ℓ ∈ Ij .

c contains points µ b1 , . . . , µ bM such that, after relabeling if necessary, Suppose M bj = µ

1 X xT,ℓ |Cj | ℓ∈C

for some nonempty Cj ⊆ Ij ,

j ∈ [M ].

Then h

bj − µj ∥2 X , E ≤ Rn,j , E ∥µ

j ∈ [M ].

The next two propositions verify the within-basin averaging hypothesis in the two merge regimes used in practice: radius merge when the number of modes is unknown, and Ward agglomerative merge when it is known. Proposition A.4 (Radius merge when the number of modes is unknown). Assume Ij ̸= ∅ for every j ∈ [M ]. Suppose the final merge is the radius-based rule with merge radius hmode . If max max ∥xT,ℓ − µj ∥ ≤

j∈[M ] ℓ∈Ij

hmode , 4

and, when the routine is applied to all k endpoints, min

5 min ∥xT,m − µj ∥ > hmode , 4

j∈[M ] m∈∪ / M j=1 Ij

b1 , . . . , µ bM such that then, for all sufficiently large n, the merged output contains points µ bj = µ

1 X xT,ℓ |Cj | ℓ∈C

for some nonempty Cj ⊆ Ij ,

j ∈ [M ].

Proposition A.5 (Ward agglomerative merge when the number of modes is known). Assume Ij ̸= ∅ for every j ∈ [M ]. Suppose the final merge is Ward agglomerative clustering with target number of clusters equal to M , run on an input set contained in {xT,ℓ : ℓ ∈ ∪M j=1 Ij }. If hmode maxj∈[M ] maxℓ∈Ij ∥xT,ℓ − µj ∥ ≤ 4 , then for all sufficiently large n, the merged output contains b1 , . . . , µ bM such that points µ bj = µ

1 X xT,ℓ |Cj | ℓ∈C

for some nonempty Cj ⊆ Ij ,

j ∈ [M ].

It remains to pass from basinwise endpoints to the merged estimator. The proof below first upgrades the local stay-in-basin event to all basin-started trajectories and then invokes Lemma 22 together with the appropriate merge proposition. Proof of Theorem 4.5. Recall from Definition 6 that M \ \

global Gn,T = Aglobal ∩ n

stay,j En,T (ℓ).

j=1 ℓ∈Ij

By Lemma 21, Pr(Aglobal ) ≥ 1 − Cglobal,stat n−2 for all sufficiently large n. We focus on Aglobal . n n global local,j Since An ⊆ An and x0,ℓ ∈ B(µj , rj ) for every ℓ ∈ Ij , the verification carried out in the proof of Theorem 4.3 shows that, under the same tuning n , log n

m≍

T = CT log n,

h = hopt ,

the hypotheses of Proposition A.3 hold for the trajectory {xt,ℓ }Tt=0 for each pair (j, ℓ). Therefore

stay,j (ℓ) X , Aglobal ≥ 1 − n−5 Pr En,T n

for all sufficiently large n. Since M X

|Ij | ≤ k ≍ log n,

j=1

a union bound gives 

Pr

M \ \

 stay,j  ≥ 1 − kn−5 . En,T (ℓ) X , Aglobal n

j=1 ℓ∈Ij

Hence global Pr(Gn,T ) ≥ 1 − Cglobal n−2

for some constant Cglobal > 0 and all sufficiently large n. Now fix j ∈ [M ] and ℓ ∈ Ij . Since Aglobal ⊆ Alocal,j , the initialization satisfies x0,ℓ ∈ B(µj , rj ), n n global and An does not involve the ascent-stage randomness, Proposition A.2 applies conditionally on global (X , An ). Thus, for some constant Cj > 0,

E ∥xT,ℓ − µj ∥ 1E stay,j (ℓ) X , Aglobal n n,T 2

2(β−1)

≤ (1 − κj η) ∥x0,ℓ − µj ∥ + Cj h

log n 1 + + + ηdσ 2 . d+2 nh mhd+2

Using the generic bound (16), the optimized bound (17), and the same transient-term choice as ⋆ , this in Theorem 4.3, with T = CT log n, m ≍ n/ log n, h = hopt , and CT ≥ CT⋆ := maxj∈[M ] CT,j becomes

E ∥xT,ℓ − µj ∥ 1E stay,j (ℓ) X , Aglobal n n,T 2

≤ CnonDP,j

log n 2(β−1) d+2β

+ CDP,j

T d polylog(n, δ) β−1 d+β

n2 ε2modes

global stay,j Since Gn,T ⊆ En,T (ℓ), we have

global = E ∥xT,ℓ − µj ∥2 X , Gn,T

E ∥xT,ℓ − µj ∥2 1G global X , Aglobal n n,T

global Gn,T

X , Aglobal n

≤

E ∥xT,ℓ − µj ∥2 1E stay,j (ℓ) X , Aglobal n

n,T

global Gn,T

X , Aglobal n

global Using Pr Gn,T X , Aglobal ≥ 1 − kn−5 , and absorbing (1 − kn−5 )−1 into the constants for all n sufficiently large n, we obtain

global ≤ CnonDP,j E ∥xT,ℓ − µj ∥2 X , Gn,T

log n 2(β−1) d+2β

+ CDP,j

T d polylog(n, δ) β−1 d+β

n2 ε2modes

for every j ∈ [M ] and every ℓ ∈ Ij . global Therefore Lemma 22 applies with E = Gn,T . If the number of modes is unknown and radius merge is used, Proposition A.4 supplies the required within-basin representation. If the number of modes is known and Ward agglomerative merge is used, Proposition A.5 supplies it. In either c contains points µ b1 , . . . , µ bM satisfying the claimed conditional MSE bound. This proves the case, M theorem.

A.4

Minimax lower bound

Proof of Theorem 4.6. Let p0 be a C ∞ density on Rd with a unique mode at 0, such that p0 (0) > 1/C and ∇2 log p0 (0) ≺ 0. Clearly p0 ∈ Hβ (L) for all β. We next create the localized shift alternatives as follows. R R For a function ψ ∈ Hβ (L) with ψ(0) = 1, ψ(u)du = 0, |ψ(u)|du ≤ C ,∥∇ψ(0)∥ ≥ 1/C , let us define pθ (x) := p0 (x)(1 + θhβ ψ(x/h)). for θ ∈ {+1, −1}. Since translations preserve Hölder regularity, for h small enough we have pθ ∈ Hβ (L). Recall that sp (x) = ∇ log p(x). For any x, using a Taylor expansion of ∇ log p0 around x, we have θhβ−1 1 ∇p0 (x) + ∇ψ(x/h) spθ (x) = pθ (x) ψ(x/h)

Since ∥∇2 log p0 (x)∥ ≤ C for sufficiently small h, it can be checked that inf

x:∥x∥≤h

∥sp1 (x) − sp−1 (x)∥ ≥ Chβ−1 .

(18)

Note that Z

|ψ(x/h)|dx = h

TV(p1 , p−1 ) = h

|ψ(u)|du ≤ Chβ+d

p1 (x) · (hβ ψ(x/h))2 dx

β+d

and Z

KL(p1 , p−1 ) =

p1 (x) p1 (x) log dx ≤ 2 p−1 (x) 2β+d

≤ Ch

p1 (uh)(ψ(u))2 du ≤ Ch2β+d .

We next use Lemma 6.1 of Karwa and Vadhan (2018) to write #

Pθk (Tε,δ ∈ S) − δ ′ ≤ ε′ sup log Pθk′ (Tε,δ ∈ S) S;k,k′ ∈{−1,1};Pθ (Tε,δ ∈S)>δ ′ k

where δ ′ = exp(ε′ )nδTV(p1 , p−1 ), ε′ = 6εnTV(p1 , p−1 ), and for a fixed x, let Tε,δ (x) be any ε, δ transcript from n iid observations. Next, following lemmas C4 to C7 of Cai et al. (2024) we have: T

ε,δ ))2 ≤ 4 min{ε′ (exp(ε′ ) − 1), nKL(p1 , p−1 )} + 32 exp(2ε′ )n2 δ 2 (TV(p1 , p−1 ))2 (TV(p1 ε,δ , p−1

≤ C min{n2 ε2 h2β+2d , nh2β+d } + n2 δ 2 h2β+2d Now choosing h = max{n

1 − 2β+d

1 − β+d

, (nε)

} and since δ = o(n−1 ) by assumption, we have T

ε,δ ))2 ≤ c (TV(p1 ε,δ , p−1

for a sufficiently small constant c > 0. Now using Le Cam two-point lemma, in particular Lemma 1 of Yu (1997), in conjunction with (18) we have for any x with ∥x∥ ≤ h that 2(β−1)

− 2β+d

inf sup E∥sb(x) − s(x)∥2 ≥ n b s

+ (nε)

−

2(β−1) β+d

Now to obtain rates for mode estimation, note that, if sb is any score estimator and estimated mode b, we have is x b) = sb(x b) − s(x b) + s(x b) = sb(x b) − s(x b) + s(0) + (H(ξ))(x b − 0) 0 = sb(x b for some t ∈ [0, 1], which implies for ξ = tx b − 0∥2 ≥ [Eλmax (H(ξ))]−2 E∥sb(x b) − s(x b)∥2 . E∥x

where H(ξ) := ∇s(x)|x=ξ is the Hessian of log density evaluated at ξ. By assumptions, the Hessian at all x has bounded singular values, and hence the minimax lower bound for score estimation implies b − x0 (P )∥2 ≥ C inf sup E∥x b x

inf

sup E∥sb(x) − sP (x)∥2 ≳ n

b s,∥x∥≤h P

2(β−1)

− 2β+d

+ (nε)

−

2(β−1) β+d

A.5

Proofs of auxiliary lemmas

This subsection verifies the supporting lemmas used earlier. The proofs are grouped to mirror the main argument: privacy first, then local geometric and analytic control, then stabilization and minibatch bounds, and finally the DAP and global lemmas. Differential Privacy lemmas. Proof of Lemma 1. Fix x ∈ Rd . Write p := pbX (x),

p′ := pbX ′ (x),

g := ∇pbX (x),

g ′ := ∇pbX ′ (x).

Since X and X ′ differ in one entry, the KDE and KDE-gradient differences satisfy |p − p′ | ≤

2K∞ , nhd

∥g − g ′ ∥ ≤

2GK . nhd+1

Now set d′ := max{p′ , pfloor }.

d := max{p, pfloor },

Then d, d′ ≥ pfloor , and since u 7→ max{u, pfloor } is 1-Lipschitz, |d − d′ | ≤ |p − p′ |. Using sbA,pfloor ;X (x) =

clipA (g) , d

sbA,pfloor ;X ′ (x) =

clipA (g ′ ) , d′

we obtain clipA (g) − clipA (g ′ ) 1 1 + clipA (g ′ ) − d d d′ 1 ∥ clipA (g ′ )∥ ≤ ∥ clipA (g) − clipA (g ′ )∥ + |d − d′ |. pfloor d d′

sbA,pfloor ;X (x) − sbA,pfloor ;X ′ (x) ≤

The clipping map clipA is the Euclidean projection onto the closed ball of radius A, hence it is 1-Lipschitz, and also ∥ clipA (g ′ )∥ ≤ A. Therefore sbA,pfloor ;X (x) − sbA,pfloor ;X ′ (x) ≤

∥g − g ′ ∥ +

pfloor Sh (A, pfloor ) = . n

A 2GK −(d+1) 2AK∞ −d + |p − p′ | ≤ h h 2 n pfloor pfloor n p2floor

This proves the claim. Proof of Lemma 2. We prove privacy of the full T -round ascent stage. Let us fix a round t ∈ {0, . . . , T − 1}, and condition on the full transcript up to the start of round t. Under this conditioning, the current iterate xt is fixed. The round-t mechanism releases Yt = sbA,pfloor ;Bt (xt ) + Zt ,

Zt ∼ N (0, σ 2 Id ).

Let Bt and Bt′ be neighboring minibatches of size m, differing in one entry. Applying Lemma 1 with n replaced by m gives sbA,pfloor ;Bt (xt ) − sbA,pfloor ;Bt′ (xt ) ≤

Sh (A, pfloor ) . m

Thus the full-data ℓ2 -sensitivity for one round is ∆t =

Sh (A, pfloor ) . m

We write ρ :=

m , n

ε∗ := log 1 +

eεiter − 1 n(eεiter − 1) = log 1 + , ρ m

δ ∗ :=

δiter nδ = . ρ 2mT

By the Gaussian mechanism (Dwork et al., 2014, Theorem 3.22), adding Gaussian noise with standard deviation r 1.25 ∆t σ = ∗ 2 log ∗ ε δ makes the corresponding full-data round mechanism (ε∗ , δ ∗ )-DP. Since 1.25 · 2mT 2.5mT 1.25 = = , ∗ δ nδ nδ this is exactly Sh (A, pfloor )/m σ= log 1 + n(eεiter − 1)/m

2 log

2.5mT , nδ

which is precisely (6). Now sample the minibatch uniformly without replacement. By privacy amplification by subsampling without replacement (Balle et al., 2018), the actual round-t mechanism is (εiter , δiter )-DP. We now compose the T rounds. By advanced composition (Dwork et al., 2010), εcomp ≤

2T log(2/δ) εiter + T εiter (eεiter − 1),

δcomp ≤

δ + T δiter . 2

Since δiter = δ/(2T ), we have δcomp ≤

δ δ +T · = δ. 2 2T

Also, by the definition of εiter ,

εmodes . 2 Further, εiter ≤ log(1 + (e − 1)m/n) ≤ 1, so eεiter − 1 ≤ 2εiter , and therefore q

2T log(2/δ) εiter ≤

T εiter (eεiter − 1) ≤ 2T ε2iter ≤ 2T · where we also used εiter ≤

εmodes /(4T ). Hence εcomp ≤ εmodes .

Thus the full T -round ascent stage is (εmodes , δ)-DP. 58

εmodes εmodes = , 4T 2

Proof of Lemma 3. Let us fix the initialization pool I = {x0,1 , . . . , x0,k }., a round t ∈ {0, . . . , T −1}, and condition on the full transcript up to the start of round t. Under this conditioning, the current iterates xt,1 , . . . , xt,k are fixed. Let us define k

Kt := C̄h (xt,ℓ , xt,r ) ℓ,r=1 . Let Bt and Bt′ be neighboring minibatches of size m, differing in one entry. Define the concurrent stabilized score difference ∆t := st (Bt ) − st (Bt′ ) ∈ Rkd . Also define the corresponding density and gradient difference vectors k

∆pt := pbBt (xt,ℓ ) − pbBt′ (xt,ℓ )

ℓ=1

∆gt := ∇pbBt (xt,ℓ ) − ∇pbBt′ (xt,ℓ )

and

k ℓ=1

By Proposition 8 of Hall et al. (2013), applied with the exponential kernel C̄h , the corresponding RKHS sensitivities satisfy 1/2

1/2

∥∆pt ∥HC̄ ≤ h

2I0 h−d , m

Let us now define Ψ(u, v) :=

∥∆gt ∥Hd ≤ C̄h

clipA (v) , max{u, pfloor }

2I1 h−(d+1) . m

u ∈ R, v ∈ Rd .

For any u, u′ ∈ R and v, v ′ ∈ Rd , ∥Ψ(u, v) − Ψ(u′ , v ′ )∥ ≤

1 pfloor

∥v − v ′ ∥ +

A |u − u′ |. 2 pfloor

Applying this componentwise across the k starts and using (a + b)2 ≤ 2a2 + 2b2 , we obtain ∥∆t ∥Hd

√ ≤ 2

C̄h

√

1 pfloor

∥∆gt ∥Hd + C̄h

1/2

A p2floor

∥∆pt ∥HC̄

1/2

I1 AI h−(d+1) + 2 0 h−d pfloor pfloor

≤

2 2 m

∆h,corr (A, pfloor ) . m

For the finite set of current evaluation points {xt,1 , . . . , xt,k }, the RKHS norm induced by C̄h corresponds to the Mahalanobis seminorm associated with the kernel matrix

Kt = C̄h (xt,ℓ , xt,r ) ℓ,r=1 ,

and

∥v∥2Hd = v ⊤ (K†t ⊗ Id )v,

v ∈ Rkd ,

C̄h

where K†t denotes the Moore–Penrose pseudoinverse. Hence the full-data round mechanism Yt = st (Bt ) + Ξt ,

Ξt ∼ N 0, σ 2 (Kt ⊗ Id ) ,

is a Gaussian mechanism calibrated to the sensitivity ∆h,corr (A, pfloor )/m measured in this seminorm. We write ρ :=

m , n

ε∗ := log 1 +

n(eεiter − 1) eεiter − 1 = log 1 + , ρ m

δ ∗ :=

δiter nδ = . ρ 2mT

By the Gaussian mechanism (Dwork et al., 2014, Theorem 3.22), choosing ∆h,corr (A, pfloor )/m σ= ε∗

1.25

2 log

δ∗

makes the full-data round mechanism (ε∗ , δ ∗ )-DP. Since 1.25/δ ∗ = 2.5mT /(nδ), this is exactly (9). By privacy amplification by subsampling without replacement (Balle et al., 2018), the actual round-t mechanism is (εiter , δiter )-DP. By advanced composition (Dwork et al., 2010), the joint mechanism (Y0 , . . . , YT −1 ) is (εmodes , δ)-DP. f := {xT,1 , . . . , xT,k } is a deterministic function of the fixed Now the pre-merge endpoint set M f is also (εmodes , δ)initialization pool I and the joint mechanism output (Y0 , . . . , YT −1 ). Hence M c is obtained from M f by deterministic merging, M c is (εmodes , δ)-DP DP. Since the final estimator M by post-processing (Proposition 2.1 of Dwork et al., 2014). A.5.1

Local geometric and analytic lemmas.

Proof of Lemma 4. Fix j ∈ [M ]. Since p ∈ Σ(β, Lj ; Uj ) with β > 2, the function p is C 2 on Uj . In particular, p, ∇p, and ∇2 p are continuous in a neighborhood of µj . By Assumption 2(i), p(µj ) > 0. Hence, by continuity of p, there exists an open neighborhood Vj ⊂ Uj of µj such that p(µj ) >0 for all x ∈ Vj . p(x) ≥ 2 Therefore log p is well defined on Vj , and on Vj we have ∇2 log p(x) =

∇2 p(x) ∇p(x)∇p(x)⊤ − . p(x) p(x)2

Since p, ∇p, and ∇2 p are continuous on Vj , and p(x) is bounded away from zero there, it follows that x 7→ ∇2 log p(x) is continuous at µj . Hence, for every ξ > 0, there exists rej (ξ) > 0 such that ∇2 log p(x) − ∇2 log p(µj ) ≤ ξ,

∀x ∈ B(µj , rej (ξ)).

Now suppose Assumption 4 also holds, so that ∇2 log p(µj ) ⪯ −αj Id . Taking ξ = αj /2, for any x ∈ B(µj , rej (αj /2)), ∇2 log p(x) ⪯ ∇2 log p(µj ) + ∇2 log p(x) − ∇2 log p(µj ) I ⪯ −αj I + This proves the claim.

αj αj I = − I. 2 2

Proof of Lemma 5. Fix j ∈ [M ] and x ∈ B(µj , rj ). Since ∇ log p(µj ) = 0, the fundamental theorem of calculus gives Z 1

∇2 log p µj + s(x − µj ) (x − µj ) ds.

∇ log p(x) =

Therefore ⟨x − µj , ∇ log p(x)⟩ =

Z 1

(x − µj )⊤ ∇2 log p µj + s(x − µj ) (x − µj ) ds.

By Assumption 4, since rj ≤ rej and hence ∇2 log p(y) ⪯ −

αj I 2

for all y ∈ B(µj , rj ).

Hence ⟨x − µj , ∇ log p(x)⟩ ≤ −

αj ∥x − µj ∥2 . 2

If x ∈ ∂B(µj , rj ), then ∥x − µj ∥ = rj , so ⟨x − µj , ∇ log p(x)⟩ ≤ −

αj rj2 . 2

Proof of Lemma 6. Fix j ∈ [M ] and write Bj := B(µj , rj ). Fix s ∈ {0, 1, 2} and a multi-index α with |α| = s. Recall n x − X 1 X i ∂ α pb(x) = ∂αK . d+s nh h i=1 For each x ∈ Bj , decompose

∂ α pb(x) − ∂ α p(x) = ∂ α pb(x) − E ∂ α pb(x) + E ∂ α pb(x) − ∂ α p(x) . We first bound the bias term. A change of variables gives E ∂ α pb(x) =

1 hs

∂ α K(u) p(x − hu) du.

Since Bj ⊂ Uj and Bj is compact, there exists ηj > 0 such that the closed ηj -neighborhood of Bj is contained in Uj . Split the integral as α

E ∂ pb(x) =

1 hs

1 ∂ K(u) p(x−hu) du+ s h ∥u∥≤ηj /h

Z ∥u∥>ηj /h

∂ α K(u) p(x−hu) du =: I1,α (x)+I2,α (x).

For the local term I1,α (x), if ∥u∥ ≤ ηj /h and x ∈ Bj , then x − hu ∈ Uj . Let ℓ = ⌊β⌋ ≥ 2. Since p ∈ Σ(β, Lj ; Uj ), Taylor’s theorem with integral remainder yields, uniformly for x ∈ Bj and ∥u∥ ≤ ηj /h, p(x − hu) =

X (−h)|γ| |γ|≤ℓ

γ!

(∂ γ p)(x) uγ + Rℓ (x, u),

|Rℓ (x, u)| ≤ Cj hβ ∥u∥β .

Inserting this into I1,α (x) gives I1,α (x) =

X (−h)|γ|−s

γ!

|γ|≤ℓ

(∂ γ p)(x)

∂ α K(u) uγ du +

∥u∥≤ηj /h

1 hs

Z ∥u∥≤ηj /h

∂ α K(u) Rℓ (x, u) du.

For each |γ| ≤ ℓ, write Z

∂ α K(u) uγ du =

∂ α K(u) uγ du −

∥u∥≤ηj /h

∂ α K(u) uγ du.

∥u∥>ηj /h

Since K is a kernel of order ℓ, the full-space moments satisfy Z

∂ α K(u) uγ du =

  (−1)|α| α!,

γ = α,

 0,

γ ̸= α, |γ| ≤ ℓ,

by integration by parts and the order-ℓ moment conditions. Hence I1,α (x) = ∂ α p(x) + R1,α (x) + R2,α (x), where R1,α (x) = −

X (−h)|γ|−s

γ!

|γ|≤ℓ

and

1 hs

R2,α (x) =

(∂ γ p)(x)

∂ α K(u) uγ du

∥u∥>ηj /h

Z ∥u∥≤ηj /h

∂ α K(u) Rℓ (x, u) du.

Since Bj is compact and p ∈ Σ(β, Lj ; Uj ), all derivatives ∂ γ p with |γ| ≤ ℓ are bounded on Bj . Using Assumption 1(iv), for each |γ| ≤ ℓ and |α| ≤ 2, Z

|∂ α K(u)| ∥u∥|γ| du ≤

h β−|γ| Z

ηj

∥u∥>ηj /h

∥u∥β |∂ α K(u)| du,

since |γ| ≤ ℓ ≤ β. Therefore |R1,α (x)| ≤ Chβ−s ,

x ∈ Bj .

Also, using the remainder bound and Assumption 1(iv), β−s

|R2,α (x)| ≤ Cj h

∥u∥β |∂ α K(u)| du ≤ Chβ−s ,

x ∈ Bj .

For the tail term I2,α (x), using Assumption 2(iii) and Assumption 1(iv), |I2,α (x)| ≤

pmax hs

|∂ α K(u)| du ≤

∥u∥>ηj /h

pmax h β hs ηj

∥u∥β |∂ α K(u)| du ≤ Chβ−s ,

uniformly in x ∈ Bj . Combining the bounds for I1,α (x) and I2,α (x) yields sup E ∂ α pb(x) − ∂ α p(x) ≤ C hβ−s .

x∈Bj

We next bound the stochastic term. Define the centered process Zα (x) =

n 1X fx (Xi ) − Efx (X) , n i=1

x − u

fx (u) := h−(d+s) ∂ α K

By Assumption 1(iii), ∥fx ∥∞ ≤ h−(d+s) sup |∂ α K(u)| < ∞. u∈Rd

Also, by Assumption 2(iii) and Assumption 1(iv), Efx (X)2 =

1 h2(d+s)

x − y 2

∂αK

p(y) dy ≤ h−(d+2s) pmax

Z Rd

(∂ α K(u))2 du ≲ h−(d+2s) ,

uniformly over x ∈ Bj , after the change of variables u = (x − y)/h. Bernstein’s inequality therefore gives, for any fixed x ∈ Bj and any t > 0, !

c nt2 Pr |Zα (x)| > t ≤ 2 exp − −(d+2s) . h + h−(d+s) t

If ∥x − x′ ∥ ≤ ρ, then by the mean value theorem and Assumption 1(iii), |fx (u) − fx′ (u)| ≤ C h−(d+s+1) ∥x − x′ ∥, uniformly in u, and hence |Zα (x) − Zα (x′ )| ≤ C h−(d+s+1) ρ. d Let {xm }N m=1 be a ρ-net of Bj with N ≲ (rj /ρ) . Then

sup |Zα (x)| ≤ max |Zα (xm )| + C h−(d+s+1) ρ. 1≤m≤N

x∈Bj

Choose

t := C

log n , nhd+2s

ρ := c hd+s+1 t,

so that the second term is at most t/2. A union bound over the net points and the above Bernstein bound yield, for all sufficiently large n, !

Pr sup |Zα (x)| > t

≤ n−5 .

x∈Bj

Consequently, with probability at least 1 − n−5 , s

sup |Zα (x)| ≲

x∈Bj

log n . nhd+2s

Combining the bias and stochastic bounds yields, with probability at least 1 − n−5 , s

sup ∂ α pb(x) − ∂ α p(x) ≲ hβ−s +

x∈Bj

log n . nhd+2s

For fixed s, there are finitely many multi-indices α with |α| = s; a union bound over them absorbs into constants. A final union bound over s ∈ {0, 1, 2} gives, for all sufficiently large n, with probability at least 1 − n−4 , 

sup ∥∇s pb(x) − ∇s p(x)∥ ≤ C hβ−s +

x∈Bj



log n  , nhd+2s

s = 0, 1, 2,

as claimed.

Proof of Lemma 7. Fix j ∈ [M ] and set Bj := B(µj , rj ). By Assumption 2(i), p(x) ≥ pmin,j > 0 for all x ∈ Bj . By Lemma 6 with s = 0, with probability at least 1 − n−4 , s

sup |pb(x) − p(x)| ≤ C hβ +

x∈Bj

log n . nhd

Under Assumption 3, h → 0 and nhd+4 / log n → ∞. Since h ≤ 1 for all sufficiently large n, these q imply hβ + log n/(nhd ) = o(1). Hence, for all sufficiently large n, the right-hand side is at most pmin,j /2. On this event, 1 1 inf pb(x) ≥ inf p(x) − sup |pb(x) − p(x)| ≥ pmin,j − pmin,j = pmin,j . x∈Bj x∈Bj 2 2 x∈Bj Thus the claim holds with cj := pmin,j /2 and probability at least 1 − n−4 for sufficiently large n. Proof of Lemma 8. Fix j ∈ [M ] and write Bj := B(µj , rj ). Let ℓ = log p and ℓb = log pb. Define the events

A :=

inf pb(x) ≥ cj ,

x∈Bj

B :=

 2  \ s=0



sup ∥∇s pb(x) − ∇s p(x)∥ ≤ C hβ−s +

x∈Bj

 log n  . nhd+2s 

Lemma 7 gives Pr(A) ≥ 1 − n−4 for all large n. Lemma 6 gives Pr(B) ≥ 1 − n−4 for all large n. Therefore, by the union bound, Pr(Ac ∪ B c ) ≤ Pr(Ac ) + Pr(B c ) ≤ 2n−4 , so Pr(A ∩ B) ≥ 1 − 2n−4 .

(19)

Work on the event A∩B below; in particular inf x∈Bj pb(x) ≥ cj and also p(x) ≥ pmin,j for all x ∈ Bj . For any x ∈ Bj , the mean value theorem gives b |ℓ(x) − ℓ(x)| = | log pb(x) − log p(x)| =

|pb(x) − p(x)| , ξ(x)

where ξ(x) lies between pb(x) and p(x). Hence ξ(x) ≥ min{cj , pmin,j } and b sup |ℓ(x) − ℓ(x)| ≤ x∈Bj

1 sup |pb(x) − p(x)|. min{cj , pmin,j } x∈Bj

On B, the right-hand side is bounded by Cj

hβ +

log n nhd

For any x ∈ Bj , ∇pb(x) ∇p(x) ∇pb(x) − ∇p(x) b ∇ℓ(x) − ∇ℓ(x) = − = + ∇p(x) b p(x) p(x) pb(x)

1 1 . − b p(x) p(x)

Since |pb−1 − p−1 | = |pbpb−p| , on A and p ≥ pmin,j , p b ∥∇ℓ(x) − ∇ℓ(x)∥ ≤

1 ∥∇p(x)∥ ∥∇pb(x) − ∇p(x)∥ + |pb(x) − p(x)|. cj cj pmin,j

Taking suprema over x ∈ Bj and using boundedness of ∥∇p(x)∥ on Bj (from Assumption 2) gives a constant Cj such that ! b − ∇ℓ(x)∥ ≤ Cj sup ∥∇ℓ(x) x∈Bj

sup ∥∇pb(x) − ∇p(x)∥ + sup |pb(x) − p(x)| . x∈Bj

x∈Bj

On B,each term is bounded by the stated rate with s = 1 and s = 0, hence the RHS is bounded by Cj hβ−1 + 2b

log n nhd+2

. For any x ∈ Bj , ∇2 pb(x) ∇2 p(x) − p(x) pb(x)

∇ ℓ(x) − ∇ ℓ(x) =

−

∇pb(x)∇pb(x)⊤ ∇p(x)∇p(x)⊤ − p(x)2 pb(x)2

For the Hessian-fraction term, ∇2 pb − ∇2 p 1 1 ∇2 pb ∇2 p − = + ∇2 p − , pb p pb pb p

so on A and p ≥ pmin,j , ∇2 pb(x) ∇2 p(x) 1 ∥∇2 p(x)∥ − ≤ ∥∇2 pb(x) − ∇2 p(x)∥ + |pb(x) − p(x)|. pb(x) p(x) cj cj pmin,j For the quadratic-gradient term, add and subtract ∇pb(x)∇p(x)⊤ : ∇pb∇pb⊤ ∇p∇p⊤ ∇pb(∇pb − ∇p)⊤ (∇pb − ∇p)∇p⊤ 1 1 − = + + ∇p∇p⊤ 2 − 2 , 2 2 2 2 pb p pb pb pb p

and use ∥uv ⊤ ∥ ≤ ∥u∥ ∥v∥ and 1 pb2

−

1 |p − pb| |p + pb| (pmax,j + ∥pb∥∞ ) = ≤ |pb − p| 2 2 2 b p p p c2j p2min,j

∇pb(x)∇pb(x)⊤ ∇p(x)∇p(x)⊤ ≤ Cj − sup pb(x)2 p(x)2 x∈Bj

sup ∥∇pb(x) − ∇p(x)∥ + sup |pb(x) − p(x)| .

x∈Bj

Combine the two displays and take suprema over x ∈ Bj to obtain ! 2b

sup ∥∇ ℓ(x)−∇ ℓ(x)∥ ≤ Cj

sup ∥∇ pb(x) − ∇ p(x)∥ + sup ∥∇pb(x) − ∇p(x)∥ + sup |pb(x) − p(x)| .

x∈Bj

On B, each supremum is bounded by the corresponding rate with s = 2, 1, 0, hence the RHS is bounded by Cj hβ−2 +

log n nhd+4

On A ∩ B we have, simultaneously for s = 0, 1, 2,  b sup ∥∇s ℓ(x) − ∇s ℓ(x)∥ ≤ Cj hβ−s + x∈Bj



log n  . nhd+2s

Together with (19) this proves the lemma.

Proof of Lemma 9. Let A be the event in Lemma 8, let B be the event in Lemma 7, and let C be the event that there exists a deterministic constant Dj > 0 such that n 1 X x − Xi K sup d nh h x∈B(µj ,rj ) i=1

n 1 X x − Xi sup ∇K d nh h x∈B(µj ,rj ) i=1

≤ Dj ,

≤ Dj .

The same net-and-Bernstein argument used in the proof of Lemma 6, applied to the bounded integrable classes x−u 2 x−u 2 , u 7→ K , u 7→ ∇K h h yields Pr(C) ≥ 1 − n−4 for all sufficiently large n. On A ∩ B ∩ C, every requirement in the definition of Alocal,j is satisfied. Therefore n A ∩ B ∩ C ⊆ Alocal,j . n Hence, for all sufficiently large n, Pr Alocal,j ≥ Pr(A ∩ B ∩ C) ≥ 1 − Pr(Ac ) − Pr(B c ) − Pr(C c ) n

≥ 1 − 2n−4 − n−4 − n−4 = 1 − 4n−4 . This proves the claim.

A.5.2

Stabilization and minibatch-control lemmas.

Proof of Lemma 10. Fix j ∈ [M ] and work on Alocal,j . For every x ∈ B(µj , rj ) and every i ∈ [n], n ∥gi (x)∥ =

hd+1

∇K

x − Xi h

≤ GK h−(d+1) .

Also, by Definition 5, pb(x) ≥ cj =

pmin,j pmin ≥ 2 2

Hence ∥qi (x)∥ =

for all x ∈ B(µj , rj ).

∥gi (x)∥ 2GK −(d+1) h = C∗ , ≤ pb(x) pmin

uniformly over x ∈ B(µj , rj ) and i ∈ [n]. Proof of Lemma 11. For fixed x ∈ B(µj , rj ), write ∆p,t (x) := pbBt (x) − pb(x),

∆g,t (x) := ∇pbBt (x) − ∇pb(x).

These are centered finite-population sample means. By the standard variance formula for sampling without replacement, h

E |∆p,t (x)|2 x, X , Alocal,j ≤ n i

≤ E ∥∆g,t (x)∥2 x, X , Alocal,j n where

1 x − Xi Ki (x) := d K , h h

n n−m 1X Ki (x)2 , · m(n − 1) n i=1 n n−m 1X ∥Gi (x)∥2 , · m(n − 1) n i=1

x − Xi Gi (x) := d+1 ∇K . h h 1

Now, n n 1X 1 1 X x − Xi Ki (x)2 = d K d n i=1 h nh i=1 h

and

2 #

n n 1 1 X x − Xi 1X 2 ∥Gi (x)∥ = d+2 ∇K d n i=1 h nh i=1 h

, 2#

By Definition 5, both bracketed quantities are bounded by Dj uniformly over x ∈ B(µj , rj ) on Alocal,j . Therefore, n h

E |∆p,t (x)|2 x, X , Alocal,j ≤ n

Dj , mhd

after enlarging constants if necessary. Next define f (u, v) :=

E ∥∆g,t (x)∥2 x, X , Alocal,j ≤ n

clipA (v) . max{u, pfloor } 67

Dj mhd+2

By Proposition A.1, on Alocal,j and for all sufficiently large n, n

f pb(x), ∇pb(x) = ∇ log pb(x)

for all x ∈ B(µj , rj ).

Moreover, exactly as in the proof of Lemma 1, ∥f (u, v) − f (u′ , v ′ )∥ ≤

1 pfloor

∥v − v ′ ∥ +

A |u − u′ |. 2 pfloor

Hence ∥sbA,pfloor ;Bt (x) − ∇ log pb(x)∥2

= f pbBt (x), ∇pbBt (x) − f pb(x), ∇pb(x) ≤

2 p2floor

2A2 |∆p,t (x)|2 . p4floor

∥∆g,t (x)∥2 +

Taking conditional expectation and using the previous bounds gives h

E ∥sbA,pfloor ;Bt (x) − ∇ log pb(x)∥2 x, X , Alocal,j ≤ n

Cj,ζ mhd+2

after enlarging constants, since h−d ≤ h−(d+2) for all sufficiently large n. Finally, Jensen’s inequality gives the mean bound.

Proof of Lemma 12. For fixed x, use the notation from the proof of Lemma 11. As in that proof, for fixed x these are centered finite-population sample means and, on Alocal,j , n n 1X Ki (x)2 ≲ h−d , n i=1

n 1X ∥Gi (x)∥2 ≲ h−d−2 , n i=1

while the deterministic envelopes are |Ki (x)| ≤ K∞ h−d ,

∥Gi (x)∥ ≤ GK h−(d+1) .

Bernstein-type concentration inequalities for sampling without replacement therefore yield constants c, C > 0 such that

Pr |∆p,t (x)| > u

x, X , Alocal,j n

cmu2 ≤ 2 exp − −d , h + h−d u

and, for each coordinate r ∈ [d],

Pr |(∆g,t (x))r | > u

x, X , Alocal,j n

Set

un := C

log(eT n) , mhd+2 68

cmu2 ≤ 2 exp − −d−2 . h + h−d−1 u

un vn := √ . d

Because mhd+2 / log(eT n) → ∞, we have hun = C{h2 log(eT n)/(mhd+2 )}1/2 = o(1), and similarly hvn = o(1). Thus the linear Bernstein terms are lower order. After increasing C if necessary,

≤ (eT n)−7 , Pr |∆p,t (x)| > un x, X , Alocal,j n and, for each coordinate r ∈ [d],

≤ d−1 (eT n)−7 Pr |(∆g,t (x))r | > vn x, X , Alocal,j n for all sufficiently large n. A union bound gives

≤ (eT n)−7 . Pr ∥∆g,t (x)∥ > un x, X , Alocal,j n Now define f (u, v) :=

clipA (v) . max{u, pfloor }

By Proposition A.1, f pb(x), ∇pb(x) = ∇ log pb(x) on Alocal,j for all sufficiently large n, and n

∥f (u, v) − f (u′ , v ′ )∥ ≤

1 A ∥v − v ′ ∥ + 2 |u − u′ |. pfloor pfloor

Hence, on the event {|∆p,t (x)| ≤ un , ∥∆g,t (x)∥ ≤ un }, we have 1

∥sbA,pfloor ;Bt (x) − ∇ log pb(x)∥ ≤

pfloor

A p2floor

un .

Absorbing constants into Cj,mb proves the first display. The second follows by conditioning at each time t and applying a union bound over t = 0, . . . , T − 1. Proof of Lemma 13. Set s

rn,1 := hβ−1 +

log n , nhd+2

rn,2 := hβ−2 +

log n . nhd+4

By Lemma 8, on Alocal,j , n sup

b ∥∇ℓ(x) − ∇ℓ(x)∥ ≤ Cj,1 rn,1 ,

sup

x∈B(µj ,rj )

b ∥∇2 ℓ(x) − ∇2 ℓ(x)∥ ≤ Cj,2 rn,2

x∈B(µj ,rj )

for deterministic constants Cj,1 , Cj,2 > 0. Fix xt ∈ B(µj , rj ), write δt := xt −µj , and recall that ∇ℓ(µj ) = 0. By the fundamental theorem of calculus, Z 1

b j) + kn (xt ) = ∇ℓ(µ

b j + sδt ) ds δt . ∇2 ℓ(µ

Adding and subtracting ∇2 ℓ, and using Lemma 5, gives ⟨δt , kn (xt )⟩ ≤ −

αj b j )∥ + Cj,2 rn,2 ∥δt ∥2 . ∥δt ∥2 + ∥δt ∥ ∥∇ℓ(µ 2 69

b j )∥ ≤ Cj,1 rn,1 . Therefore Since ∇ℓ(µj ) = 0, Lemma 8 also gives ∥∇ℓ(µ

⟨δt , kn (xt )⟩ ≤ −

αj ∥δt ∥2 + Cj,1 rn,1 ∥δt ∥ + Cj,2 rn,2 ∥δt ∥2 . 2

Applying Young’s inequality to the middle term and using 2 rn,1 ≲ h2(β−1) +

log n nhd+2

yields the first display. For the second display, let Hj := supx∈B(µj ,rj ) ∥∇2 ℓ(x)∥. The same expansion gives Z 1

∇ ℓ(µj + sδt ) ds δt + Cj,2 rn,2 ∥δt ∥ ≤ Cj,1 rn,1 + Hj ∥δt ∥ + Cj,2 rn,2 ∥δt ∥.

b j )∥ + ∥kn (xt )∥ ≤ ∥∇ℓ(µ 0

Squaring both sides and using 2 rn,1 ≲ h2(β−1) +

log n , nhd+2

2 rn,2 ≲ h2(β−2) +

log n , nhd+4

gives the second display.

Proof of Lemma 14. Because Bt is a uniformly random m-subset of {1, . . . , n}, Pr(i ∈ Bt ) =

m , n

Pr(i, j ∈ Bt ) =

m(m − 1) n(n − 1)

(i ̸= j).

Writing the minibatch mean with indicators, n 1 X qi (xt )1{i ∈ Bt }, q̄t = m i=1

we obtain E[q̄t | xt , X , Alocal,j ]= n

n 1X qi (xt ) = ∇ log pb(xt ) =: kn (xt ). n i=1

For the second moment, n 1 mX m(m − 1) X 2 ∥q (x )∥ + qi (xt )⊤ qj (xt ) . i t m2 n i=1 n(n − 1) i̸=j

E[∥q̄t ∥

| xt , X , Alocal,j ]= n

Using X i̸=j

qi⊤ qj =

n X

i=1

−

n X

∥qi ∥2 = n2 ∥kn (xt )∥2 −

i=1

n X

∥qi (xt )∥2 ,

i=1

we get E[∥q̄t ∥2 | xt , X , Alocal,j ]= n

n n−m X n−m ∥qi (xt )∥2 + 1 − ∥kn (xt )∥2 . mn(n − 1) i=1 m(n − 1)

(20)

On Alocal,j , Definition 5 gives pb(xt ) ≥ cj > 0, and there exists a constant Dj > 0 such that n n x − Xi 1 X ∇K d h x∈B(µj ,rj ) nh i=1

≤ Dj .

sup

Since qi (xt ) =

gi (xt ) , pb(xt )

gi (xt ) :=

x − X

∇K hd+1

it follows that n n n x − X 2 x − X 2 1X 1 1 1 1X 1 X t i t i = ∇K ∥qi (xt )∥2 = ∇K · d+1 d 2 2 d+2 n i=1 h nh i=1 h pb(xt ) n i=1 h pb(xt ) h

≤

Dj −d−2 h . c2j

Thus E[∥q̄t ∥2 | xt , X , Alocal,j ]≲ n

n − m −d−2 h + ∥kn (xt )∥2 . m(n − 1)

Combining this with (20) and Lemma 13 yields the claim. Proof of Lemma 15. Set kn (xt ) = ∇ log pb(xt ). Then

rt = sbA,pfloor ;Bt (xt ) − kn (xt ) − q̄t − kn (xt ) . Hence, by (a + b)2 ≤ 2a2 + 2b2 , E[∥rt ∥2 | xt , X , Alocal,j ] n h

≤ 2E ∥sbA,pfloor ;Bt (xt ) − kn (xt )∥2 xt , X , Alocal,j + 2E[∥q̄t − kn (xt )∥2 | xt , X , Alocal,j ]. n n The first term is O((mhd+2 )−1 ) by Lemma 11. For the second term, (20) gives 2

E[∥q̄t − kn (xt )∥

| xt , X , Alocal,j ]= n

n n−m X 1 ∥qi (xt )∥2 ≲ . mn(n − 1) i=1 mhd+2

This proves the claim. Proof of Lemma 16. We condition on X and on Alocal,j throughout. By construction, the process n (j) T −1 et )t=0 is adapted to the algorithmic filtration and satisfies (x (j)

et x

∈ B(µj , rj )

(j)

for all t = 0, . . . , T − 1.

Define (j)

et ) − ∇ log p(x et ), bt := ∇ log pb(x

(j)

et ) − ∇ log pb(x et ). ζt := sbA,pfloor ;Bt (x

Then (j)

(j)

et ) − ∇ log p(x et ) + zt = bt + ζt + zt . sbA,pfloor ;Bt (x

(j)

et Since x

∈ B(µj , rj ) for all t, Definition 5 gives 

max ∥bt ∥ ≤ Cj hβ−1 +

0≤t≤T −1



log n  . nhd+2

(j)

−1 et )Tt=0 Next, Lemma 12 applies to the adapted process (x , yielding



Pr max ∥ζt ∥ ≤ Cj,mb 0≤t≤T −1



log(eT n)  ≥ 1 − T (eT n)−6 . X , Alocal,j n mhd+2

For the Gaussian noise, recall that each zt is a d-dimensional Gaussian row with covariance 2 σ Id . Standard Gaussian norm concentration and a union bound over t = 0, . . . , T − 1 give a deterministic constant C2,j > 0 such that

max ∥zt ∥ ≤ C2,j σ d log(eT n) X , Alocal,j n

0≤t≤T −1

≥ 1 − T (eT n)−6 .

By (8), √ ∆h,corr (A, pfloor ) = 2 2

1/2

I1 AI h−(d+1) + 2 0 h−d . pfloor pfloor

Since A and pfloor are fixed in n and h → 0 under Assumption 3, one has h ≤ 1 for all sufficiently large n, hence h−d ≤ h−(d+1) . Therefore ∆h,corr (A, pfloor ) ≲ h−(d+1) ≍ C∗ . √ Using (9) and the fact that εiter ≲ εmodes / T , there exists a constant C3,j > 0 such that p

C∗ T polylog(T, n, δ) σ ≤ C3,j nεmodes for all sufficiently large n. Therefore C2,j σ

C∗ d log(eT n) ≤ Ce2,j

T d polylog(T, n, δ) nεmodes

after enlarging constants. Intersecting the minibatch and Gaussian-noise good events and combining the three bounds above yields (j) max ∥bt + ζt + zt ∥ ≤ Ξn,T,m,h 0≤t≤T −1

with conditional probability at least 1 − 2T (eT n)−6 . This proves the lemma.

A.5.3

Initialization lemmas.

Proof of Lemma 17. The design gives hDAP ≍ n−1/(d+2β) (log n)1/(d+2β)+1/d → 0, ρinit

ρinit ≍ (log n)−1/d → 0.

Since minj rj > 0 and mini̸=j ∥µi − µj ∥ ≥ c0 > 0, choose n0 such that, for all n ≥ n0 , √ d ρinit 1 2 hDAP ≤ , ρinit ≤ min rj , ρinit ≤ c0 . 2 8 2 j 3 For each j, the lattice-grid property gives zr ∈ Zn with √ d ρinit ρinit hDAP ≤ < . ∥zr − µj ∥ ≤ 2 8 4 Thus zr ∈ Gj,n , proving (i). If x, a ∈ B(µj , ρinit /4), then ∥x − a∥ ≤ ∥x − µj ∥ + ∥a − µj ∥ ≤

ρinit ρinit + < ρinit , 4 4

which proves (ii). If x ∈ B(µi , ρinit /4), a ∈ B(µj , ρinit /4), and i ̸= j, then ∥x − a∥ ≥ ∥µi − µj ∥ − ∥x − µi ∥ − ∥a − µj ∥ ≥ c0 −

ρinit ≥ ρinit , 2

which proves (iii).

Proof of Lemma 18. Fix j ∈ [M ], and write ℓ = log p. By the local Hessian assumption, there are constants 0 < aj < Aj < ∞ such that, for every 0 < t ≤ rj , inf

∥x−µj ∥≤t

ℓ(x) ≥ ℓ(µj ) − aj t2 ,

sup

ℓ(y) ≤ ℓ(µj ) − Aj t2 .

t≤∥y−µj ∥≤rj

Choose t = ρinit /4. Since p is bounded below on B(µj , rj ), the preceding display implies, after decreasing Aj − aj if necessary, that for all x ∈ B(µj , ρinit /4) and all y ∈ B(µj , rj ) \ B(µj , ρinit /4), p(x) − p(y) ≥ cj ρ2init for a constant cj > 0. For zr ∈ B(µj , ρinit /4), the change of variables x = zr + hDAP u gives mr = hdDAP

Z ∥u∥≤1

p(zr + hDAP u) du.

The same identity holds for zs . Because p ∈ C 2 (Uj ) and hDAP /ρinit → 0, sup

|p(z + hDAP u) − p(z)| ≤ ChDAP = o(ρ2init ).

z∈B(µj ,rj ) ∥u∥≤1

Hence, for zr ∈ Zn ∩ B(µj , ρinit /4) and zs ∈ Zn ∩ (B(µj , rj ) \ B(µj , ρinit /4)), mr − ms = hdDAP

≥ hdDAP

∥u∥≤1

{p(zr + hDAP u) − p(zs + hDAP u)} du {cj ρ2init − o(ρ2init )} du

≥ ∆j hdDAP ρ2init for some ∆j > 0 and all n ≥ nj . Taking the infimum over r and the supremum over s proves the claim.

Proof of Lemma 19. For each r, ur is the average of Bernoulli variables with mean mr . Since pmax < ∞, mr ≤ pmax Vol{x : ∥x − zr ∥ ≤ hDAP } ≤ ChdDAP . Bernstein’s inequality gives, for every t > 0, nt2 Pr(|ur − mr | > t) ≤ 2 exp − 2mr + 2t/3

≤ 2 exp −

nt2 ChdDAP + t

Taking t = C1

s  hd



log Ncand  DAP log Ncand +  n n



with C1 large enough and union bounding over r = 1, . . . , Ncand yields

max

1≤r≤Ncand

−4 |ur − mr | > t ≤ CNcand .

−4 d 2 Since Ncand ≤ Ch−d DAP ≤ Cn, the right side is at most Cn . It remains to compare t with hDAP ρinit . Under the DAP design,

hdDAP ≍ Therefore

log n n

hdDAP log Ncand /n hdDAP ρ2init

and

d/(d+2β)

ρ2init ≍ (log n)−2/d .

≍ n−β/(d+2β) (log n)1/2−d/(2(d+2β))+2/d → 0

log Ncand /n ≍ n−2β/(d+2β) (log n)1−d/(d+2β)+2/d → 0. d 2 hDAP ρinit

Thus, for any fixed cconc > 0, t ≤ cconc hdDAP ρ2init for all large enough n. This proves the stated bound.

Proof of Lemma 20. Fix m ∈ [M ]. Let zr ∈ Cm,n , and suppose that [

zr ∈ /

B(µj , rj ).

j∈[M ]\{m}: p(µj )≥p(µm )

Since zr ∈ / B(µm , rm ) by definition of Cm,n , the point zr lies outside the basin neighborhood of every mode whose height is at least p(µm ). By the compactness separation argument used for the DAP competitive regions, there is a constant ηm > 0 such that, outside these higher-or-equal modal neighborhoods, p(zr ) ≤ p(µm ) − ηm . On the other hand, for every zs ∈ Gm,n , the definition of Gm,n gives ∥zs − µm ∥ ≤ ρinit /4. Since ρinit → 0 and p is continuous at µm , there is nm such that, for all n ≥ nm , p(zs ) ≥ p(µm ) − ηm /4

for every zs ∈ Gm,n .

Using the local-mass expansion from Lemma 18, there is a constant Cm > 0 such that, uniformly over the grid points under consideration, ms − hdDAP Vol(B(0, 1))p(zs ) ≤ Cm hd+1 DAP . Increase nm if needed so that

ηm Vol(B(0, 1)). 4

2Cm hDAP ≤ Then, for every zs ∈ Gm,n ,

ms − mr ≥ hdDAP Vol(B(0, 1)){p(zs ) − p(zr )} − 2Cm hd+1 DAP η m d ≥ hDAP Vol(B(0, 1)) p(µm ) − − (p(µm ) − ηm ) − 2Cm hd+1 DAP 4 ηm ≥ Vol(B(0, 1))hdDAP . 2 Increase nm if needed so that γ0 ρ2init <

ηm Vol(B(0, 1)). 2

Then mr <

inf

s: zs ∈Gm,n

ms − γ0 hdDAP ρ2init ,

contradicting zr ∈ Cm,n . Therefore Cm,n ⊆

[ j∈[M ]\{m}: p(µj )≥p(µm )

B(µj , rj ).

It remains to cover this set. Since Cm,n ⊆ Zn ⊆ Q, it is enough to cover the fixed public box √ Q. A bounded box in Rd can be covered by at most CQ ρ−d init cubes of side length ρinit /(2 d), and each such cube is contained in a Euclidean ball of radius ρinit /2. Hence N (Cm,n , ∥ · ∥, ρinit /2) ≤ CQ ρ−d init . Since ρinit ≍ (log n)−1/d and M ≥ 1, there is a constant C ′ > 0 such that ′ CQ ρ−d init ≤ C M log n.

Absorbing C ′ into Lcap proves the covering bound. A.5.4

Global-control lemmas.

Proof of Lemma 21. For each j ∈ [M ], let Aj := Alocal,j , and let B := Einit . By Lemma 9, for each n j ∈ [M ], Pr(Acj ) ≤ 4n−4 for all sufficiently large n. By Proposition 4.4, there exists Cinit,cov > 0 such that Pr(B c ) ≤ Cinit,cov n−2 for all sufficiently large n. Therefore, for all sufficiently large n, Pr

M \

Aj ∩ B ≥ 1 −

j=1

M X

Pr(Acj ) − Pr(B c ) ≥ 1 − 4M n−4 − Cinit,cov n−2 .

j=1

Hence Pr Aglobal ≥ 1 − Cglobal,stat n−2 for a constant Cglobal,stat > 0 and sufficiently large n. n

f to M: c once a merged point The next argument isolates the only point needed to pass from M is an average of endpoints attached to a single population mode, the pre-merge rate is inherited unchanged.

Proof of Lemma 22. By convexity of ∥ · ∥2 , 2

1 X b j − µj ∥ = (xT,ℓ − µj ) ∥µ |Cj | ℓ∈C 2

≤

1 X ∥xT,ℓ − µj ∥2 . |Cj | ℓ∈C j

Taking conditional expectation given (X , E) and using the assumed endpoint bounds for all ℓ ∈ Cj ⊆ Ij gives h i i 1 X h b j − µ j ∥2 X , E ≤ E ∥µ E ∥xT,ℓ − µj ∥2 X , E ≤ Rn,j . |Cj | ℓ∈C j

Proof of Proposition A.4. Assume first that ℓ, m ∈ Ij . Then ∥xT,ℓ − xT,m ∥ ≤ ∥xT,ℓ − µj ∥ + ∥xT,m − µj ∥ ≤

hmode < hmode . 2

So all endpoints from the same basin lie within the merge radius of one another. If ℓ ∈ Ii , m ∈ Ij , and i ̸= j, then ∥xT,ℓ − xT,m ∥ ≥ ∥µi − µj ∥ − ∥xT,ℓ − µi ∥ − ∥xT,m − µj ∥ ≥ ∥µi − µj ∥ −

hmode . 2

Since ∆min := mini̸=j ∥µi − µj ∥ > 0 is fixed and hmode → 0, we have ∆min −

hmode > hmode 2

for all sufficiently large n. Thus endpoints from different basins are not merged. If the routine is applied to all k endpoints, the additional assumption min

5 min ∥xT,m − µj ∥ > hmode 4

j∈[M ] m∈∪ / M j=1 Ij

implies that no endpoint outside ∪M j=1 Ij can lie within distance hmode of any endpoint in Ij . Therefore, for each j ∈ [M ], radius merge produces a cluster consisting only of indices from Ij ; call that index set Cj . Its merged point is exactly bj = µ

1 X xT,ℓ , |Cj | ℓ∈C j

with nonempty Cj ⊆ Ij . This proves the proposition. Proof of Proposition A.5. Let Nmax := maxj∈[M ] |Ij |. If A, B ⊆ Ij are nonempty subclusters from the same mode, then their centroids satisfy ∥x̄A − x̄B ∥ ≤

hmode , 2

so their Ward linkage obeys ∆Ward (A, B) =

|A||B| Nmax h2mode ∥x̄A − x̄B ∥2 ≤ . |A| + |B| 4

If A ⊆ Ii and B ⊆ Ij with i ̸= j, then ∥x̄A − x̄B ∥ ≥ ∥µi − µj ∥ −

hmode hmode ≥ ∆min − , 2 2

|A||B| where ∆min := mini̸=j ∥µi − µj ∥ > 0. Since |A|+|B| ≥ 12 for nonempty clusters,

∆Ward (A, B) ≥

1 hmode 2 ∆min − . 2 2 77

Because Nmax ≤ k ≍ log n and hmode → 0, one has Nmax h2mode → 0. Hence, for all sufficiently large n, hmode 2 1 Nmax h2mode ∆min − . < 4 2 2 So every within-mode Ward merge has smaller cost than every cross-mode Ward merge. Since the algorithm is stopped at M clusters, the final partition is mode-pure. If Cj denotes the final cluster associated with µj , then Cj ⊆ Ij , and the corresponding Ward centroid is bj = µ

1 X xT,ℓ . |Cj | ℓ∈C j

This proves the proposition.

Additional Experimental Results

The code to reproduce all experimental results can be found at https://github.com/ArkaB-DS/ DP-GRAMS. All experiments are implemented in Python and run on a machine with an ARM CPU (8 cores, 8 logical processors), 8.6 GB RAM, and macOS 15.6.1. Unless otherwise noted, reported summaries are averages over 20 runs with standard errors, and runtimes are reported in seconds.

B.1

Private Mode Estimation

This section provides additional diagnostics for the two mode-estimation benchmarks used in the paper: the bivariate 4-modal Gaussian mixture introduced in Section 3.1 and the bivariate 5-modal t-mixture studied in Section 5.3. These results supplement the main-text privacy–utility curves in Figures 1 and 3 by showing run-to-run variability, sensitivity to the clipping threshold C∗ , minibatch size m, and step size η, and MSE and runtime summaries across (n, ε). Bivariate 4-modal Gaussian mixture. These diagnostics assess whether the four-corners Gaussian results in Section 3.1 are stable across independent private runs and moderate tuning changes. Figures 11–13 show that, in the Gaussian benchmark, the private estimates remain concentrated near the four modal basins and do not exhibit sharp degradation near the selected clipping, minibatch, or step-size defaults. Table 2 shows that most of the MSE reduction occurs when moving from the tightest privacy budget to moderate ε, with additional improvement as n increases. The reported runtimes are small in this implementation because DP-GRAMS follows a fixed number of private starts rather than running mean shift from every data point. These diagnostics support the main privacy–utility trends shown in Figure 1. Bivariate 5-modal t-mixture. We repeat the same diagnostics for the heavier-tailed, heterogeneous mixture studied in Section 5.3.

Mode Estimation Across Runs for 4-Modal Bivariate Gaussian Mixture

Run 1

Run 2

Run 3

Run 4

Run 5

Run 6

Run 7

Run 8

Run 9

Run 10

Run 11

Run 12

Run 13

Run 14

Run 15

Run 16

Run 17

Run 18

Run 19

Run 20

True modes

DP-GRAMS

Figure 11: Bivariate 4-modal Gaussian mixture: grid of 20 DP-GRAMS runs on one fixed dataset. Each subplot shows KDE contours with true modes (green stars), non-private mean-shift estimates (blue crosses), and DP-GRAMS estimates (orange circles). The grid visualizes run-to-run variability from private DAP initialization and injected ascent noise.

MSE vs m across n for 4-Modal Bivariate Gaussian Mixture

MSE vs C * across n for 4-Modal Bivariate Gaussian Mixture

0.5

n = 700

n=700 n=1000 n=2000 n=5000

0.4 MSE

n = 1000

1.0 MSE

0.6

0.5 0.0 10

0.3

n = 2000

0.2

n = 5000

1.0 MSE

0.1 0.0 0.1

0.2

0.3

0.4

0.5

Clip Multiplier (C * scaling)

0.6

0.7

0.5 0.0

0.8

Minibatch size m

(b) MSE vs. minibatch size m.

(a) MSE vs. clipping multiplier clip multiplier.

Figure 12: Hyperparameter sensitivity for DP-GRAMS on the 4-modal Gaussian mixture. (a) Effect of clipping multiplier clip multiplier on MSE for n ∈ {700, 1000, 2000, 5000} at fixed ε = 1. (b) Effect of minibatch size m on MSE across the same sample sizes. The MSE is relatively stable around the default choices for both tuning parameters. Figures 14–16 and Table 3 show the same qualitative pattern as the main t-mixture experiment: errors are large at the tightest privacy budgets, decrease sharply as ε increases, and continue to improve with n. The gap to the mean-shift baseline remains larger than in the Gaussian benchmark, reflecting the heavier tails and heterogeneous component scales. Across both synthetic benchmarks, the appendix diagnostics support the main-text privacy–utility trends and show that the selected clipping, minibatch, and step-size defaults are not isolated tuning choices.

B.2

Private Modal Regression

This subsection provides additional diagnostics for the private modal-regression experiments in Section 5.4. We first report a complementary three-component piecewise-constant design, then give expanded privacy–utility tables and hyperparameter sweeps for both this design and the sinusoidal two-component mixture from the main text. Unless otherwise stated, privacy is calibrated as in Section 5.4, and regression error is the oracle modal MSE in (15). Runtime summaries for PMS and DP-PMS are aggregated over the same 20 runs. 3-component mixture. This example is a three-component mixture with piecewise-constant conditional modes. For sample size n, let n1 = ⌊n/3⌋,

n2 = ⌊n/3⌋,

n3 = n − n1 − n2 .

We generate independent samples X1,i ∼ Uniform(0, 0.5),

Y1,i ∼ N (3, σ 2 ), 80

i = 1, . . . , n1 ,

MSE vs across n for 4-Modal Bivariate Gaussian Mixture n = 700

0.6

n = 1000

MSE

0.4 0.2 0.0 10

n = 2000

0.6

n = 5000

MSE

0.4 0.2 0.0 10

Step size

Figure 13: Step-size sensitivity for DP-GRAMS on the 4-modal Gaussian mixture. The figure reports MSE versus step size η across n ∈ {700, 1000, 2000, 5000} at fixed ε = 1. The sweep does not show sharp degradation near the selected default. X2,i ∼ Uniform(0.4, 0.7),

Y2,i ∼ N (2, σ 2 ),

i = 1, . . . , n2 ,

X3,i ∼ Uniform(0.6, 1),

Y3,i ∼ N (1, σ 2 ),

i = 1, . . . , n3 ,

with σ = 0.2. Because the predictor intervals overlap, the population conditional mode set is piecewise: 

M(x) =

  {3},       {3, 2},   

{2},

     {2, 1},      {1},

x < 0.4,

0.4 ≤ x ≤ 0.5, 0.5 < x < 0.6, 0.6 ≤ x ≤ 0.7, x > 0.7.

For this support-limited design, we use a binned conditional-DAP initializer. This is a variant of Algorithm 4: instead of selecting sparse predictor locations directly from the full public design, we partition the fixed predictor domain into bins and, within each nonempty bin, sample a private Y -anchor from a public Y -grid using an exponential-mechanism score based on nearby responses. The resulting (x, y) anchors are passed to DP-PMS as initialization points. After private ascent, a fixed-radius cleanup is applied using only the already private DP-PMS outputs, so this cleanup is post-processing. 81

Table 2: Bivariate 4-modal Gaussian mixture: MSE and runtime (mean ± SE) for DP-GRAMS and non-private mean shift across (n, ε).

MSE

Runtime (s)

DP-GRAMS

700

0.1 0.25 0.5 1.0 5.0

0.7375 ± 0.0978 0.2166 ± 0.0233 0.0903 ± 0.0201 0.0590 ± 0.0070 0.0551 ± 0.0072

0.00851 ± 0.00000 0.00851 ± 0.00000 0.00851 ± 0.00000 0.00851 ± 0.00000 0.00851 ± 0.00000

0.00250 ± 0.00002 0.00253 ± 0.00001 0.00251 ± 0.00002 0.00248 ± 0.00001 0.00256 ± 0.00003

0.09364 ± 0.00014 0.09364 ± 0.00014 0.09364 ± 0.00014 0.09364 ± 0.00014 0.09364 ± 0.00014

1000

0.1 0.25 0.5 1.0 5.0

0.4724 ± 0.0567 0.1232 ± 0.0171 0.0596 ± 0.0087 0.0463 ± 0.0070 0.0439 ± 0.0046

0.00492 ± 0.00000 0.00492 ± 0.00000 0.00492 ± 0.00000 0.00492 ± 0.00000 0.00492 ± 0.00000

0.00331 ± 0.00005 0.00306 ± 0.00000 0.00308 ± 0.00003 0.00292 ± 0.00001 0.00296 ± 0.00002

0.17251 ± 0.00050 0.17251 ± 0.00050 0.17251 ± 0.00050 0.17251 ± 0.00050 0.17251 ± 0.00050

2000

0.1 0.25 0.5 1.0 5.0

0.1887 ± 0.0210 0.0539 ± 0.0062 0.0307 ± 0.0028 0.0222 ± 0.0019 0.0207 ± 0.0023

0.00421 ± 0.00000 0.00421 ± 0.00000 0.00421 ± 0.00000 0.00421 ± 0.00000 0.00421 ± 0.00000

0.00482 ± 0.00006 0.00473 ± 0.00001 0.00472 ± 0.00001 0.00470 ± 0.00000 0.00473 ± 0.00001

0.69787 ± 0.00188 0.69787 ± 0.00188 0.69787 ± 0.00188 0.69787 ± 0.00188 0.69787 ± 0.00188

5000

0.1 0.25 0.5 1.0 5.0

0.0588 ± 0.0062 0.0184 ± 0.0023 0.0121 ± 0.0010 0.0107 ± 0.0015 0.0090 ± 0.0015

0.00194 ± 0.00000 0.00194 ± 0.00000 0.00194 ± 0.00000 0.00194 ± 0.00000 0.00194 ± 0.00000

0.01137 ± 0.00022 0.01107 ± 0.00001 0.01144 ± 0.00027 0.01110 ± 0.00001 0.01105 ± 0.00001

6.09457 ± 0.12074 6.09457 ± 0.12074 6.09457 ± 0.12074 6.09457 ± 0.12074 6.09457 ± 0.12074

Figure 17 shows that PMS separates the piecewise modal branches, whereas LOWESS averages across them and misses the conditional multimodality. At ε = 1, DP-PMS tracks the branch structure after private response-direction ascent. Table 4 shows a steep privacy–utility transition: for example, at n = 2000, DP-MSE drops from 0.6380 at ε = 0.1 to 0.000174 at ε = 1. Figure 18 shows that performance does not sharply degrade near the selected clipping and minibatch defaults. These diagnostics indicate that the binned conditional-DAP initializer is useful for this supportlimited, piecewise-constant modal-regression design. Sinusoidal 2-mixture. We extend the sinusoidal two-component experiment in Section 5.4 by reporting the full privacy–utility table and the clipping and minibatch sensitivity sweeps. The datagenerating mechanism, privacy grid, and oracle loss are the same as in the main-text experiment. 82

Mode Estimation Across Runs for 5-Modal Bivariate t-Mixture

Run 1

Run 2

Run 3

Run 4

Run 5

Run 6

Run 7

Run 8

Run 9

Run 10

Run 11

Run 12

Run 13

Run 14

Run 15

Run 16

Run 17

Run 18

Run 19

Run 20

True modes

DP-GRAMS

Figure 14: Bivariate 5-modal t-mixture: grid of 20 DP-GRAMS runs on one fixed dataset. Each subplot shows KDE contours with true modes (green stars), non-private mean-shift estimates (blue crosses), and DP-GRAMS estimates (orange circles). The grid visualizes run-to-run variability from private DAP initialization and injected ascent noise in the heavier-tailed setting.

MSE vs m across n for 5-Modal Bivariate t-Mixture

MSE vs C * across n for 5-Modal Bivariate t-Mixture

n = 700

MSE

0.5 0.4 MSE

n = 1000

1.5

n=700 n=1000 n=2000 n=5000

1.0 0.5 0.0

0.3

n = 2000

0.2

n = 5000

1.5 MSE

0.1

1.0 0.5

0.0 0.4

0.5

0.6

0.7

0.8

Clip Multiplier (C * scaling)

0.9

0.0

1.0

Minibatch size m

(b) MSE vs. minibatch size m.

(a) MSE vs. clipping multiplier clip multiplier.

Figure 15: Hyperparameter sensitivity for DP-GRAMS on the 5-modal t-mixture. (a) Effect of clipping multiplier clip multiplier on MSE for n ∈ {700, 1000, 2000, 5000} at fixed ε = 1. (b) Effect of minibatch size m on MSE across the same sample sizes. The sweeps do not show sharp degradation near the selected defaults, despite the heavier tails and heterogeneous component scales. Table 5 summarizes oracle DP-MSE, PMS-MSE, and runtime across the sinusoidal privacy– utility grid. The table supports the main trend in Figure 4: the largest DP-MSE reductions occur between the smallest privacy budgets and moderate ε, while for larger n further privacy-budget increases yield smaller improvements and the private error moves closer to the PMS-MSE scale. Figure 19 shows that the sinusoidal experiment is not sharply sensitive to moderate changes in the clipping multiplier or minibatch size near the selected defaults. Together, the piecewise-constant and sinusoidal designs show that DP-PMS can recover both support-limited branches and smooth nonlinear modal curves under the fixed-design responseprivacy setup. The two examples differ in difficulty: the three-component design shows a sharper drop in DP-MSE at larger n and ε, while the sinusoidal design retains a visible gap from PMS-MSE even at the largest privacy budgets. In both cases, the clipping and minibatch sweeps indicate that the reported results are not driven by a narrow tuning choice.

B.3

Private Clustering on Simulated Data

This subsection complements the blobs clustering experiment in Section 5.5.1 by providing additional numerical summaries and hyperparameter sweeps for DP-GRAMS-C. Figure 6 reports privacy–utility curves in ARI, NMI, and centroid MSE versus ε across n ∈ {700, 1000, 2000, 5000}. Table 6 shows that DP-GRAMS-C improves sharply from ε = 0.1 to moderate privacy budgets and then largely stabilizes in ARI and NMI, while centroid MSE remains small once ε is moderate. DP-k-Means also improves with ε, but is generally weaker in ARI and NMI except at the loosest 84

MSE vs across n for 5-Modal Bivariate t-Mixture n = 700

n = 1000

0.8 MSE

0.6 0.4 0.2 0.0 10

n = 2000

n = 5000

0.8 MSE

0.6 0.4 0.2 0.0 10

Step size

Modal Regression for 3-Component Mixture Data PMS LOWESS

0.0

0.2

0.4

0.6

0.8

1.0

4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0

Privacy-Utility Tradeoff for 3-Component DP Modal Regression

Differentially Private Modal Regression for 3-Component Mixture Data DP-PMS ( =1)

MSE

4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0

Figure 16: Step-size sensitivity for DP-GRAMS on the 5-modal t-mixture. The figure reports MSE versus step size η across n ∈ {700, 1000, 2000, 5000} at fixed ε = 1. The sweep does not show sharp degradation near the selected default, although the heavier-tailed setting exhibits some variability across sample sizes.

0.0

0.2

0.4

0.6

0.8

1.0

PMS n=200 PMS n=500 PMS n=1000 PMS n=2000 DP-PMS n=200 DP-PMS n=500 DP-PMS n=1000 DP-PMS n=2000

Privacy budget

Figure 17: Private modal regression on three-component piecewise-constant mixture data. Panels (a) and (b) use one representative dataset with n = 500; panel (b) uses (ε, δ) = (1, 10−5 ). (a) PMS tracks the modal branches, whereas LOWESS smooths across them. (b) DP-PMS recovers the piecewise branch structure under privacy. (c) Privacy–utility tradeoff: oracle MSE in (15) versus ε on a log scale for n ∈ {200, 500, 1000, 2000}.

privacy budgets and largest sample sizes. We then fix ε = 1 and study sensitivity to the clipping threshold C∗ and minibatch size m. Effects of minibatch size m and clipping threshold C∗ . Figure 20 shows how clustering

Table 3: Bivariate 5-modal t-mixture: MSE and runtime (mean ± SE) for DP-GRAMS and nonprivate mean shift across (n, ε).

MSE

Runtime (s)

DP-GRAMS

700

0.1 0.25 0.5 1.0 5.0

2.8695 ± 0.2815 1.4499 ± 0.2419 0.6244 ± 0.0723 0.2226 ± 0.0191 0.0502 ± 0.0053

0.00450 ± 0.00000 0.00450 ± 0.00000 0.00450 ± 0.00000 0.00450 ± 0.00000 0.00450 ± 0.00000

0.00539 ± 0.00002 0.00533 ± 0.00000 0.00533 ± 0.00000 0.00532 ± 0.00000 0.00535 ± 0.00001

0.09284 ± 0.00010 0.09284 ± 0.00010 0.09284 ± 0.00010 0.09284 ± 0.00010 0.09284 ± 0.00010

1000

0.1 0.25 0.5 1.0 5.0

3.0136 ± 0.8413 1.0509 ± 0.0954 0.3243 ± 0.0381 0.1272 ± 0.0133 0.0530 ± 0.0046

0.01001 ± 0.00000 0.01001 ± 0.00000 0.01001 ± 0.00000 0.01001 ± 0.00000 0.01001 ± 0.00000

0.00654 ± 0.00002 0.00649 ± 0.00001 0.00647 ± 0.00001 0.00648 ± 0.00001 0.00659 ± 0.00004

0.16933 ± 0.00033 0.16933 ± 0.00033 0.16933 ± 0.00033 0.16933 ± 0.00033 0.16933 ± 0.00033

2000

0.1 0.25 0.5 1.0 5.0

1.1756 ± 0.1011 0.3725 ± 0.0363 0.1076 ± 0.0119 0.0405 ± 0.0020 0.0180 ± 0.0019

0.00331 ± 0.00000 0.00331 ± 0.00000 0.00331 ± 0.00000 0.00331 ± 0.00000 0.00331 ± 0.00000

0.01094 ± 0.00008 0.01083 ± 0.00001 0.01079 ± 0.00001 0.01079 ± 0.00001 0.01080 ± 0.00001

0.68341 ± 0.00030 0.68341 ± 0.00030 0.68341 ± 0.00030 0.68341 ± 0.00030 0.68341 ± 0.00030

5000

0.1 0.25 0.5 1.0 5.0

0.6145 ± 0.0600 0.1063 ± 0.0082 0.0310 ± 0.0032 0.0215 ± 0.0024 0.0125 ± 0.0015

0.00293 ± 0.00000 0.00293 ± 0.00000 0.00293 ± 0.00000 0.00293 ± 0.00000 0.00293 ± 0.00000

0.02399 ± 0.00017 0.02373 ± 0.00002 0.02372 ± 0.00001 0.02371 ± 0.00001 0.02378 ± 0.00004

4.49957 ± 0.06095 4.49957 ± 0.06095 4.49957 ± 0.06095 4.49957 ± 0.06095 4.49957 ± 0.06095

quality varies with C∗ across the sample-size grid. Figure 21 shows how performance varies with minibatch size across the same sample-size grid.

B.4

Private Clustering on Real Datasets

This subsection supplements the real-data clustering experiments in Section 5.5. We report additional Digits results, together with MNIST privacy–utility summaries and sensitivity analyses. Digits uses the default public DAP grid in a six-dimensional PCA representation, whereas MNIST uses public auxiliary candidates in a whitened five-dimensional PCA representation to avoid the high-dimensional DAP-grid bottleneck. The following subsection reports the corresponding Cancer RNA-Seq diagnostics.

MSE vs C * across n for 3-Component Modal Regression

MSE vs m across n for 3-Component DP Modal Regression

MSE

n = 200

10 n=200 n=500 n=1000 n=2000

Clip Multiplier (C * scaling)

MSE

n = 500

n = 1000

(a) MSE vs. clipping multiplier clip multiplier across n.

n = 2000

Minibatch size m

(b) MSE vs. minibatch size m across n.

Figure 18: Hyperparameter sensitivity for DP-PMS on the three-component mixture. (a) Effect of the clipping multiplier clip multiplier on oracle MSE for n ∈ {200, 500, 1000, 2000} at fixed ε = 1. (b) Effect of minibatch size m on oracle MSE across the same sample sizes. The sweeps do not show sharp degradation near the selected defaults. MSE vs C * across n for Sinusoidal Modal Regression 1

n = 200

MSE

2 × 10

MSE vs m across n for Sinusoidal 2-Component DP Modal Regression n=200 n=600 n=1200 n=2400

2 × 10

6 × 10

n = 600

n = 1200

MSE

6 × 10

Clip Multiplier (C * scaling)

2 × 10

6 × 10

(a) MSE vs. clipping multiplier clip multiplier across n.

Minibatch size m

n = 2400

Minibatch size m

(b) MSE vs. minibatch size m across n.

Figure 19: Hyperparameter sensitivity for DP-PMS on sinusoidal two-component mixture data. (a) Effect of the clipping multiplier clip multiplier on oracle MSE for n ∈ {200, 600, 1200, 2400} at fixed ε = 1. (b) Effect of minibatch size m on oracle MSE across the same sample sizes. The sweeps do not show sharp degradation near the selected defaults. Figure 22 shows that DP-GRAMS-C preserves much of the cluster structure seen under nonprivate mean shift in the PCA visualization. Figure 23 and Table 7 show a stable advantage over DP-k-Means on this dataset: DP-GRAMS-C has substantially higher ARI and NMI and lower 87

Table 4: Private modal regression on the three-component mixture: oracle DP-MSE, PMS-MSE, and runtime summaries (mean ± SE) for DP-PMS and PMS over the same 20 runs.

DP-MSE

PMS-MSE

PMS-runtime

DP-runtime

200 200 200 200

0.1 0.2 0.5 1.0

1.3397 ± 0.2489 0.9565 ± 0.1265 0.7605 ± 0.0666 0.3653 ± 0.0547

0.0020 ± 0.0005 0.0030 ± 0.0005 0.0026 ± 0.0006 0.0028 ± 0.0007

0.0125 ± 0.0007 0.0131 ± 0.0009 0.0123 ± 0.0006 0.0126 ± 0.0008

0.0129 ± 0.0007 0.0133 ± 0.0007 0.0127 ± 0.0008 0.0137 ± 0.0010

500 500 500 500

0.1 0.2 0.5 1.0

0.7823 ± 0.0798 0.8034 ± 0.0558 0.4760 ± 0.0541 0.1402 ± 0.0351

0.0017 ± 0.0004 0.0021 ± 0.0004 0.0013 ± 0.0003 0.0015 ± 0.0004

0.0409 ± 0.0011 0.0423 ± 0.0014 0.0412 ± 0.0022 0.0419 ± 0.0014

0.0144 ± 0.0009 0.0135 ± 0.0006 0.0150 ± 0.0012 0.0166 ± 0.0011

1000 1000 1000 1000

0.1 0.2 0.5 1.0

0.7579 ± 0.0648 0.7279 ± 0.0578 0.1743 ± 0.0335 0.0066 ± 0.0033

0.0016 ± 0.0002 0.0013 ± 0.0003 0.0015 ± 0.0002 0.0016 ± 0.0003

0.1070 ± 0.0039 0.1059 ± 0.0040 0.1214 ± 0.0037 0.1053 ± 0.0026

0.0166 ± 0.0010 0.0183 ± 0.0013 0.0193 ± 0.0011 0.0170 ± 0.0008

2000 2000 2000 2000

0.1 0.2 0.5 1.0

0.6380 ± 0.0570 0.2213 ± 0.0356 0.0130 ± 0.0076 0.000174 ± 0.000026

0.0016 ± 0.0002 0.0015 ± 0.0002 0.0012 ± 0.0002 0.0014 ± 0.0002

0.3441 ± 0.0098 0.4284 ± 0.0163 0.3860 ± 0.0127 0.3767 ± 0.0085

0.0225 ± 0.0009 0.0261 ± 0.0019 0.0270 ± 0.0020 0.0246 ± 0.0012

Blobs: ARI vs C * across n n = 700

Blobs: NMI vs C * across n n = 1000

n = 700

n = 1000

NMI

0.6

0.5 0.5

0.4 10

n = 2000

n = 5000

0.5 10

Clip Multiplier (C * )

(a) ARI vs. C∗ .

Clip Multiplier (C * )

(b) NMI vs. C∗ .

n = 2000 Centroid MSE

0.6

0.5 0.4

0.5 0.0

n = 1000

1.0

n = 5000

0.7

0.6

NMI

ARI

n = 2000

0.7

n = 700 Centroid MSE

0.7

ARI

Blobs: MSE vs C * across n

n = 5000

1.0 0.5 0.0 10

Clip Multiplier (C * )

Figure 20: Effect of clipping on DP-GRAMS-C for blobs. Subsampling-effect study for C∗ ∈ {0.01, 0.1, 0.5, 1.0, 2.0} at ε = 1 across n ∈ {700, 1000, 2000, 5000}. ARI, NMI, and centroid MSE vary smoothly with C∗ , with a broad range of clipping multipliers (including the default C∗ = 1) yielding near-optimal performance. centroid MSE throughout the privacy grid. The gains are not strictly monotone in every metric,

Table 5: Private modal regression on sinusoidal two-component mixture data: oracle DP-MSE, PMS-MSE, and runtime summaries (mean ± SE) for DP-PMS and PMS over the same 20 runs.

DP-MSE

PMS-MSE

PMS-runtime

DP-runtime

200 200 200 200

0.1 0.2 0.5 1.0

3.5963 ± 0.2478 0.5476 ± 0.0423 0.1388 ± 0.0069 0.1102 ± 0.0061

0.0313 ± 0.0014 0.0305 ± 0.0017 0.0326 ± 0.0013 0.0339 ± 0.0018

0.0128 ± 0.0008 0.0146 ± 0.0010 0.0127 ± 0.0007 0.0136 ± 0.0009

0.0324 ± 0.0018 0.0394 ± 0.0015 0.0358 ± 0.0013 0.0379 ± 0.0019

600 600 600 600

0.1 0.2 0.5 1.0

0.1642 ± 0.0062 0.1020 ± 0.0057 0.0642 ± 0.0042 0.0551 ± 0.0033

0.0203 ± 0.0003 0.0213 ± 0.0005 0.0218 ± 0.0006 0.0209 ± 0.0005

0.0683 ± 0.0041 0.0677 ± 0.0036 0.0578 ± 0.0025 0.0568 ± 0.0023

0.0633 ± 0.0073 0.0507 ± 0.0029 0.0494 ± 0.0032 0.0484 ± 0.0021

1200 1200 1200 1200

0.1 0.2 0.5 1.0

0.0975 ± 0.0046 0.0571 ± 0.0041 0.0469 ± 0.0028 0.0438 ± 0.0018

0.0186 ± 0.0002 0.0179 ± 0.0002 0.0187 ± 0.0002 0.0183 ± 0.0002

0.1629 ± 0.0028 0.1569 ± 0.0033 0.1560 ± 0.0035 0.1572 ± 0.0040

0.0640 ± 0.0028 0.0658 ± 0.0027 0.0629 ± 0.0023 0.0642 ± 0.0031

2400 2400 2400 2400

0.1 0.2 0.5 1.0

0.0453 ± 0.0025 0.0429 ± 0.0017 0.0426 ± 0.0026 0.0417 ± 0.0019

0.0165 ± 0.0002 0.0165 ± 0.0001 0.0165 ± 0.0002 0.0161 ± 0.0002

0.4593 ± 0.0128 0.4567 ± 0.0119 0.4713 ± 0.0132 0.4747 ± 0.0128

0.0758 ± 0.0029 0.0792 ± 0.0030 0.0738 ± 0.0025 0.0800 ± 0.0029

Blobs: ARI vs m across n n = 700

Blobs: NMI vs m across n n = 1000

n = 700

0.70

Centroid MSE

NMI

0.70 0.68

n = 2000

n = 5000

n = 2000

NMI

0.70 2

Minibatch size m

0.70

Minibatch size m

(a) ARI vs. minibatch size m.

Minibatch size m

0.05 10

n = 2000

0.68 10

0.10

n = 5000

0.72

0.75

n = 1000

0.15

Centroid MSE

ARI

n = 700

0.72

0.75

ARI

Blobs: MSE vs m across n n = 1000

Minibatch size m

(b) NMI vs. minibatch size m.

n = 5000

0.15 0.10 0.05 10

Minibatch size m

Figure 21: Effect of minibatch size on DP-GRAMS-C for blobs. Sample-size grid for m ∈ {0.01n, 0.05n, 0.1n, 0.2n, n} at ε = 1 and C∗ = 1 across n ∈ {700, 1000, 2000, 5000}. The metrics do not show sharp deterioration near the selected minibatch default.

Table 6: Blobs (simulated): privacy–utility summary for DP-GRAMS-C and DP-k-Means across (n, ε). We report ARI, NMI, centroid MSE, and runtime as mean ± SE. n

DP-k-Means

DP-GRAMS-C

ε ARI

NMI

MSE

Time (s)

ARI

NMI

MSE

Time (s)

0.1 0.2 0.5 1.0 5.0

0.383 ± 0.025 0.569 ± 0.022 0.744 ± 0.005 0.759 ± 0.003 0.758 ± 0.002

0.476 ± 0.018 2.047 ± 0.254 0.597 ± 0.014 0.605 ± 0.121 0.707 ± 0.003 0.0985 ± 0.0136 0.715 ± 0.002 0.0589 ± 0.00590 0.715 ± 0.002 0.0500 ± 0.00450

0.00317 ± 0.00018 0.00252 ± 0.00001 0.00250 ± 0.00001 0.00259 ± 0.00002 0.00359 ± 0.00023

0.376 ± 0.027 0.377 ± 0.018 0.421 ± 0.014 0.500 ± 0.022 0.722 ± 0.008

0.466 ± 0.020 0.465 ± 0.013 0.506 ± 0.010 0.545 ± 0.016 0.693 ± 0.005

1.699 ± 0.152 1.557 ± 0.145 0.981 ± 0.0553 0.771 ± 0.0845 0.126 ± 0.0180

0.00722 ± 0.00002 0.00704 ± 0.00003 0.00719 ± 0.00008 0.00703 ± 0.00006 0.01230 ± 0.00009

0.1 0.2 1000 0.5 1.0 5.0

0.486 ± 0.027 0.643 ± 0.018 0.711 ± 0.005 0.729 ± 0.003 0.740 ± 0.002

0.539 ± 0.018 1.197 ± 0.204 0.647 ± 0.012 0.248 ± 0.0406 0.689 ± 0.004 0.0791 ± 0.00747 0.701 ± 0.002 0.0567 ± 0.00701 0.707 ± 0.002 0.0426 ± 0.00318

0.00276 ± 0.00001 0.00281 ± 0.00002 0.00275 ± 0.00001 0.00283 ± 0.00003 0.00290 ± 0.00001

0.387 ± 0.025 0.360 ± 0.021 0.466 ± 0.019 0.521 ± 0.025 0.689 ± 0.013

0.477 ± 0.017 0.454 ± 0.015 0.528 ± 0.013 0.566 ± 0.017 0.675 ± 0.008

1.877 ± 0.190 1.877 ± 0.168 1.034 ± 0.107 0.616 ± 0.0861 0.167 ± 0.0390

0.00774 ± 0.00011 0.00729 ± 0.00003 0.00703 ± 0.00003 0.00719 ± 0.00006 0.01883 ± 0.00011

0.1 0.2 2000 0.5 1.0 5.0

0.593 ± 0.021 0.713 ± 0.005 0.728 ± 0.002 0.733 ± 0.003 0.739 ± 0.002

0.606 ± 0.014 0.682 ± 0.003 0.693 ± 0.002 0.695 ± 0.002 0.700 ± 0.002

0.348 ± 0.0589 0.0923 ± 0.00986 0.0611 ± 0.00472 0.0539 ± 0.00495 0.0479 ± 0.00639

0.00445 ± 0.00012 0.00428 ± 0.00002 0.00446 ± 0.00004 0.00410 ± 0.00003 0.00410 ± 0.00001

0.378 ± 0.014 0.362 ± 0.025 0.508 ± 0.025 0.554 ± 0.023 0.715 ± 0.015

0.457 ± 0.011 1.582 ± 0.145 0.448 ± 0.016 1.394 ± 0.109 0.553 ± 0.017 0.669 ± 0.0900 0.585 ± 0.014 0.439 ± 0.0746 0.686 ± 0.009 0.0827 ± 0.0419

0.00801 ± 0.00005 0.00815 ± 0.00008 0.00796 ± 0.00006 0.00767 ± 0.00005 0.02321 ± 0.00005

0.1 0.2 5000 0.5 1.0 5.0

0.717 ± 0.007 0.740 ± 0.002 0.743 ± 0.002 0.738 ± 0.002 0.731 ± 0.003

0.687 ± 0.005 0.702 ± 0.001 0.702 ± 0.001 0.701 ± 0.001 0.697 ± 0.002

0.0907 ± 0.0132 0.0342 ± 0.00325 0.0308 ± 0.00382 0.0456 ± 0.00413 0.0682 ± 0.00548

0.00744 ± 0.00004 0.00740 ± 0.00002 0.00735 ± 0.00001 0.00747 ± 0.00002 0.00736 ± 0.00001

0.405 ± 0.024 0.495 ± 0.024 0.554 ± 0.024 0.639 ± 0.024 0.729 ± 0.013

0.483 ± 0.016 1.186 ± 0.123 0.546 ± 0.016 0.763 ± 0.0914 0.588 ± 0.015 0.494 ± 0.0791 0.638 ± 0.016 0.255 ± 0.0615 0.694 ± 0.008 0.0370 ± 0.0161

0.00962 ± 0.00004 0.00973 ± 0.00004 0.01361 ± 0.00006 0.02519 ± 0.00008 0.02862 ± 0.00007

700

Clustering Comparison on Digits Dataset MS Clustering

DP-GRAMS-C

KMeans

PC2

PC1

20 Cluster 0 Cluster 1

PC1 Cluster 2 Cluster 3

20 Cluster 4 Cluster 5

30 20

Cluster Assignments & Centroids Cluster 6 Cluster 7

DP-KMeans

PC2

Cluster 8 Cluster 9

PC1

True means Estimated modes

Figure 22: Digits dataset. Two-dimensional visualization of the six-dimensional PCA clustering representation, comparing non-private mean shift, DP-GRAMS-C, k-means, and DP-k-Means, with private methods run at ε = 1. True class centroids and estimated centroids are overlaid.

but centroid MSE is lowest at the largest privacy budget and the ARI and NMI values remain consistently separated from the DP-k-Means baseline. Hyperparameter sweeps on Digits. For Digits, we sweep C∗ ∈ {0.1, 0.5, 1, 5, 10} and minibatch fractions m/n ∈ {0.1, 0.2, 0.5, 1.0} at fixed (ε, δ) = (1, 10−5 ). For each configuration, we record ARI, NMI, centroid MSE, and runtime across repeated runs. Figures 24–25 show no sharp deterioration near the selected default settings.

Privacy-Utility: ARI (Digits)

Privacy-Utility: NMI (Digits)

0.6

0.2

0.0

0.0 10

DP-GRAMS-C DP-KMeans

0.4

MSE

NMI

ARI

DP-GRAMS-C DP-KMeans

20000

0.6

0.4

Privacy-Utility: MSE (Digits)

25000

15000 10000 5000 0

modes

Figure 23: Privacy–utility on Digits: ARI, NMI, and centroid MSE versus ε on a log scale for DP-GRAMSC and DP-k-Means, with ε ∈ {0.25, 0.5, 1, 2.5, 5}. Points show averages over 20 runs with standard-error bars.

Table 7: Digits: privacy–utility summary for DP-GRAMS-C and DP-k-Means across ε. Reported metrics are ARI, NMI, centroid MSE, and runtime, summarized as mean ± SE over 20 runs.

Algorithm

ARI

NMI

Centroid MSE

Runtime (s)

0.25 0.25

DP-GRAMS-C DP-k-Means

0.384 ± 0.009 0.036 ± 0.010

0.535 ± 0.007 0.134 ± 0.021

28.93 ± 2.47 253.96 ± 12.88

0.1465 ± 0.0022 0.0335 ± 0.0009

0.5 0.5

DP-GRAMS-C DP-k-Means

0.374 ± 0.011 0.029 ± 0.007

0.527 ± 0.009 0.118 ± 0.019

31.60 ± 3.51 260.87 ± 14.63

0.1476 ± 0.0020 0.0336 ± 0.0009

1.0 1.0

DP-GRAMS-C DP-k-Means

0.376 ± 0.011 0.046 ± 0.010

0.532 ± 0.008 0.153 ± 0.018

27.69 ± 2.11 217.67 ± 11.49

0.1401 ± 0.0009 0.0321 ± 0.0011

2.5 2.5

DP-GRAMS-C DP-k-Means

0.392 ± 0.009 0.115 ± 0.016

0.548 ± 0.008 0.262 ± 0.023

28.58 ± 2.03 204.23 ± 12.83

0.1395 ± 0.0008 0.0314 ± 0.0009

5.0 5.0

DP-GRAMS-C DP-k-Means

0.399 ± 0.010 0.136 ± 0.017

0.549 ± 0.007 0.285 ± 0.020

20.96 ± 1.86 179.29 ± 12.14

0.1408 ± 0.0012 0.0294 ± 0.0007

DP-GRAMS-C on Digits: MSE vs C *

DP-GRAMS-C on Digits: ARI vs C * 0.70

0.60

15000 10000 5000 0 0.00

0.02

0.04

0.06

Clip Multiplier (C * )

0.08

0.10

NMI

0.55

ARI

Centroid MSE

DP-GRAMS-C on Digits: NMI vs C *

0.65

20000

0.50

0.65

0.45

0.60

0.40

0.55

0.35

0.00

0.02

0.04

0.06

Clip Multiplier (C * )

0.08

0.10

0.00

0.02

0.04

0.06

Clip Multiplier (C * )

0.08

0.10

Figure 24: Digits, DP-GRAMS-C: centroid MSE, ARI, and NMI versus clipping multiplier clip multiplier at (ε, δ) = (1, 10−5 ). Points show averages over 20 runs with standard-error bars.

MNIST public-candidate sensitivity. For MNIST, the main experiment uses a stratified public auxiliary candidate set of 1000 images, projects both the public candidates and the private ex91

DP-GRAMS-C on Digits: ARI vs m

120

0.62

100

0.60

80 60

0.72

0.58 0.56

40 10

Minibatch size m

0.70 0.68

0.54

DP-GRAMS-C on Digits: NMI vs m

0.74

NMI

0.64

ARI

Centroid MSE

DP-GRAMS-C on Digits: MSE vs m 140

Minibatch size m

Figure 25: Digits, DP-GRAMS-C: centroid MSE, ARI, and NMI versus minibatch size m on a log scale at (ε, δ) = (1, 10−5 ). Points show averages over 20 runs with standard-error bars.

perimental data into a whitened five-dimensional PCA representation, and runs clustering in that shared representation. No full DAP lattice grid is built for MNIST. We sweep the clipping multiplier and minibatch size at fixed (ε, δ) = (1, 10−5 ), using the same public candidate construction as in Section 5.5.2. Table 8 summarizes the MNIST privacy–utility statistics for DP-GRAMS-C and DP-k-Means.

0.426

0.312

0.424

0.310

0.422

0.308

0.420

0.306

0.418 2

Clip Multiplier (C * )

DP-GRAMS-C on MNIST: MSE vs C *

115

Centroid MSE

0.314

DP-GRAMS-C on MNIST: NMI vs C *

0.428

NMI

ARI

DP-GRAMS-C on MNIST: ARI vs C *

110 105 100

Clip Multiplier (C * )

Figure 26: MNIST, DP-GRAMS-C: ARI, NMI, and centroid MSE versus clipping multiplier clip multiplier at (ε, δ) = (1, 10−5 ). Points show averages over 20 runs with standard-error bars.

DP-GRAMS-C on MNIST: ARI vs minibatch size

0.4300

0.318

DP-GRAMS-C on MNIST: MSE vs minibatch size 110

0.4250

0.312

0.4225

NMI

0.314 0.310

0.4200

0.308

0.4175

0.306

Centroid MSE

0.4275

0.316

ARI

DP-GRAMS-C on MNIST: NMI vs minibatch size

0.4150 10

Minibatch size m

105 100 95 90

Minibatch size m

Figure 27: MNIST, DP-GRAMS-C: ARI, NMI, and centroid MSE versus minibatch size m on a log scale at (ε, δ) = (1, 10−5 ). Points show averages over 20 runs with standard-error bars.

Together, the Digits and MNIST results provide two complementary real-image checks for the private prototype-release clustering procedure. Digits uses the default public DAP grid in a six-dimensional PCA representation and shows a clear separation between DP-GRAMS-C and 92

Table 8: MNIST: privacy–utility summary for DP-GRAMS-C and DP-k-Means across ε. Reported metrics are ARI, NMI, centroid MSE, and runtime, summarized as mean ± SE over 20 runs.

Algorithm

ARI

NMI

Centroid MSE

Runtime (s)

0.05 0.05

DP-GRAMS-C DP-k-Means

0.205 ± 0.010 0.148 ± 0.008

0.339 ± 0.007 0.277 ± 0.008

160.12 ± 7.06 348.03 ± 16.49

0.5064 ± 0.0555 0.1450 ± 0.0009

0.1 0.1

DP-GRAMS-C DP-k-Means

0.272 ± 0.007 0.202 ± 0.009

0.389 ± 0.006 0.326 ± 0.007

107.56 ± 5.08 285.23 ± 15.39

0.4765 ± 0.0335 0.1450 ± 0.0007

0.2 0.2

DP-GRAMS-C DP-k-Means

0.303 ± 0.003 0.220 ± 0.009

0.417 ± 0.003 0.346 ± 0.007

99.54 ± 4.22 219.99 ± 17.62

0.5078 ± 0.0483 0.1466 ± 0.0017

0.5 0.5

DP-GRAMS-C DP-k-Means

0.309 ± 0.003 0.246 ± 0.006

0.422 ± 0.003 0.369 ± 0.004

104.70 ± 5.53 181.13 ± 14.82

0.5120 ± 0.0417 0.1978 ± 0.0014

1.0 1.0

DP-GRAMS-C DP-k-Means

0.311 ± 0.003 0.279 ± 0.003

0.424 ± 0.003 0.395 ± 0.002

108.92 ± 5.23 148.61 ± 14.47

0.3755 ± 0.0015 0.3750 ± 0.0017

DP-k-Means across the privacy grid. MNIST uses public auxiliary candidates in a whitened fivedimensional PCA representation; there the gap is smaller, but DP-GRAMS-C improves from ε = 0.05 to moderate privacy budgets and remains competitive at ε = 1. These diagnostics support the main clustering findings in Section 5.5 and indicate that the conclusions are not driven by a narrow choice of clipping threshold or minibatch size.

B.5

Private Clustering on Cancer Gene Expression (RNA-Seq)

This subsection supplements the Cancer RNA-Seq clustering experiment in Section 5.5.3. The data are standardized gene-expression profiles, clustered in a whitened six-dimensional PCA representation. DP-GRAMS-C uses the default public DAP grid in this reduced space, with no data points supplied as DAP candidates; the public DAP box is chosen using the robust R90 rule from the implementation. Centroid MSE is computed after inverse-PCA back-projection of estimated centers to standardized gene space. Table 9 summarizes privacy–utility behavior across ε ∈ {0.5, 1, 2, 5, 10} for DP-GRAMS-C and a DP-k-Means baseline. The improvement occurs between ε = 0.5 and ε = 1, where DPGRAMS-C moves from moderate clustering quality to ARI and NMI near the non-private range; for larger ε, ARI and NMI mostly stabilize while centroid MSE continues to decrease. DP-k-Means improves with ε but remains evidently worse in ARI, NMI, and centroid MSE throughout the grid. Hyperparameter sensitivity: We evaluate sensitivity of DP-GRAMS-C to the clipping multiplier clip multiplier and minibatch size m at fixed (ε, δ) = (1, 10−5 ). We sweep C∗ in the grid {0.01, 0.02, 0.05, 0.1} and minibatch fractions m/n ∈ {0.05, 0.1, 0.2, 0.5, 1.0}, averaging each 93

Table 9: Cancer RNA-Seq after gene-wise standardization and projection to a whitened six-dimensional PCA representation: privacy–utility summary for DP-GRAMS-C and DP-k-Means across ε. Reported metrics are ARI, NMI, centroid MSE, and runtime, summarized as mean ± SE over 20 runs. ε

Algorithm

ARI

NMI

Centroid MSE

Runtime (s)

0.5

DP-GRAMS-C

0.572 ± 0.028

0.668 ± 0.020

3401.42 ± 422.59

2.445 ± 0.105

0.5

DP-k-Means

0.151 ± 0.026

0.270 ± 0.027

13813.81 ± 774.31

0.0200 ± 0.0003

1.0

DP-GRAMS-C

0.773 ± 0.017

0.839 ± 0.009

1115.66 ± 173.85

2.158 ± 0.068

1.0

DP-k-Means

0.139 ± 0.033

0.246 ± 0.031

14138.50 ± 854.17

0.0196 ± 0.0002

2.0

DP-GRAMS-C

0.794 ± 0.013

0.846 ± 0.008

1146.84 ± 277.41

2.430 ± 0.103

2.0

DP-k-Means

0.260 ± 0.029

0.366 ± 0.025

10613.41 ± 543.80

0.0193 ± 0.0002

5.0

DP-GRAMS-C

0.808 ± 0.009

0.859 ± 0.004

827.88 ± 245.25

2.644 ± 0.077

5.0

DP-k-Means

0.430 ± 0.027

0.520 ± 0.022

6986.08 ± 536.68

0.0193 ± 0.0001

10.0

DP-GRAMS-C

0.811 ± 0.009

0.859 ± 0.004

800.60 ± 247.91

2.335 ± 0.062

10.0

DP-k-Means

0.458 ± 0.027

0.560 ± 0.023

5721.49 ± 534.16

0.0188 ± 0.0001

configuration over 20 runs. Figures 28 and 29 show no sharp degradation near the selected defaults. DP-GRAMS-C on Gene: Centroid vs C *

0.825

DP-GRAMS-C on Gene: NMI vs C * 0.86

0.800

0.84

1500 1000 0.02

0.04

0.06

Clip Multiplier (C * )

0.08

0.10

NMI

0.775

2000

ARI

Centroid MSE

2500

DP-GRAMS-C on Gene: ARI vs C *

0.750

0.82

0.725

0.80

0.700

0.78

0.675

0.02

0.04

0.06

Clip Multiplier (C * )

0.08

0.10

0.02

0.04

0.06

Clip Multiplier (C * )

0.08

0.10

Figure 28: Cancer RNA-Seq, DP-GRAMS-C: centroid MSE, ARI, and NMI versus clipping multiplier clip multiplier at (ε, δ) = (1, 10−5 ). Points show averages over 20 runs with standard-error bars.

DP-GRAMS-C on Gene: Centroid vs m

DP-GRAMS-C on Gene: ARI vs m

DP-GRAMS-C on Gene: NMI vs m 0.86

0.82

5000

0.80

0.85

4000

0.78

0.84

NMI

6000

ARI

Centroid MSE

7000

3000

0.76

0.83

2000

0.74

0.82

1000

0.72 10

Minibatch size m

0.81 10

Minibatch size m

Figure 29: Cancer RNA-Seq, DP-GRAMS-C: centroid MSE, ARI, and NMI versus minibatch size m on a log scale at (ε, δ) = (1, 10−5 ). Points show averages over 20 runs with standard-error bars.

Together with the main Cancer RNA-Seq results in Figures 9 and 10, these appendix diagnostics show that DP-GRAMS-C remains effective in a high-dimensional gene-expression task after PCA reduction. The method reaches ARI and NMI near the non-private range by moderate privacy budgets, while centroid MSE continues to improve at larger ε. The sensitivity plots indicate that this behavior is not driven by a narrow choice of clipping threshold or minibatch size.

Algorithms

Algorithm 3: DP-GRAMS-C: DP-GRAMS based Clustering Input : Private data S = {Xi }ni=1 ⊂ Rd ; privacy parameters (ε, δ); optional number of clusters kest ; initialization privacy fraction p0 ; bandwidth multiplier; clipping multiplier; optional minibatch size m; optional public candidate set Zpub . Output: Private cluster centers M; optional deterministic assignments y = (y1 , . . . , yn ) obtained from the released centers. 1 Choose bandwidth h from the sample size n, dimension d, and bandwidth multiplier; 2 Run DP-GRAMS on S with parameters (ε, δ, p0 , h, m), using Zpub if supplied and otherwise constructing the public DAP grid internally, to obtain private candidate modes raw Mraw = {µraw 1 , . . . , µr }; 3 if kest is specified then 4 Merge Mraw into at most kest clusters using agglomerative clustering, and replace each cluster by its mean to obtain M; else 6 Merge nearby points in Mraw using the default distance-threshold rule to obtain M;

for i = 1, . . . , n do 8 Assign each point deterministically to its nearest released private center yi = arg minj ∥Xi − Mj ∥2 ;

return M, and y when assignments are requested for post-processing or evaluation;

Algorithm 4: DP-PMS: Differentially Private Partial Mean Shift Input : Data {(Xi , Yi )}ni=1 , with predictor locations treated as fixed and public; predictor evaluation grid X = {x(1) , . . . , x(G) }; privacy parameters (ε, δ); minibatch size m; number of iterations T ; ascent bandwidth h; DAP bandwidth hscore ; clipping multiplier cclip ; initialization privacy fraction p0 ; sparse-start multiplier κinit ; stepsize η; public response candidate grid Ypub . c (g) )}G . Output: Private conditional mode estimates {M(x g=1 1 Set εinit = p0 ε, εasc = (1 − p0 )ε, εdraw = εinit /k, clipping threshold C = cclip /h ; 2 Select predictor locations x1 , . . . , xk from the predictor design, with k = ⌈κinit log n⌉; 3 for j = 1, . . . , k do P 4 Compute the public local predictor count Nj = ni=1 1{|Xi − xj | ≤ hscore }; 5 For each response candidate z ∈ Ypub , compute the conditional local-mass utility Pn

uj (z) = 6

(0)

Sample yj

i=1 1{|Xi − xj | ≤ hscore , |Yi − z| ≤ hscore }

max{Nj , 1}

∈ Ypub using weights proportional to exp{(εdraw Nj /2)uj (z)}, z ∈ Ypub ;

7 Calibrate the Gaussian noise scale σ for T ascent steps using privacy budget (εasc , δ); 8 for t = 0, . . . , T − 1 do 9 10 11 12

Sample one minibatch Bt ⊂ [n] uniformly without replacement, with |Bt | = m; for j = 1, . . . , k do (t)

Compute localized scalar contributions  (t) (t) qi (xj , yj ) = 1 

13 14

(t)

For each i ∈ Bt , compute wij = exp −{(Xi − xj )2 + (Yi − yj )2 }/(2h2 ) ;



(t) (t) wij (Yi − yj ) (t)  wℓj ̸= 0 , P (t) ℓ∈Bt ℓ∈Bt wℓj

(t)

Form the clipped update direction ∆j =

(t)

Kt (j, ℓ) = exp−

16 17

(t) (t) i∈Bt max{−C, min{qi (xj , yj ), C}};

Form the correlation matrix Kt ∈ Rk×k over current joint states (xj , yj ): 

i ∈ Bt ,

(t)

(xj , yj ) − (xℓ , yℓ ) h

 2

Draw correlated Gaussian noise ξ (t) ∼ N (0, σ 2 Kt ); for j = 1, . . . , k do (t+1)

Update only the response coordinate: yj

(t)

= yj + η ∆j + ξj

(T ) 18 For each x(g) ∈ X , collect terminal values yj whose fixed predictor locations xj lie near

c (g) ); x(g) , and merge nearby response values by post-processing to obtain M(x c (g) )}G ; 19 return {M(x g=1

Record · ID 422232 · SHA-256 dffa74ae0a81d70f
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.