ConceptioArchivearXiv CS
arXiv CSopen access

Differential Privacy of Gaussian Process Posterior Sampling

Unknown · 2026 · arxiv_cs
arXiv CS · Papers · License: Open Access · 2026
Open Source ↗Direct PDF ↓
cryptography, security, privacy, cybersecurity

Differential Privacy of Gaussian Process Posterior Sampling

Tomasz Maciazek School of Mathematics, University of Bristol

arXiv:2606.17995v1 [stat.ML] 16 Jun 2026

Abstract We study the privacy of releasing posterior sample paths from a Gaussian process (GP) when the entire training set including covariates and responses is private. Unlike standard differential-privacy (DP) mechanisms that add external noise, posterior sampling is random by construction. We show that this intrinsic randomness yields DP guarantees by deriving explicit Rényi-DP bounds for GP posterior sample-path release. The bounds separate posterior-mean leakage from data-dependent posterior-covariance leakage showing that meaningful privacy depends sharply on effective ridge regularisation. We apply membership-inference attacks to show that empirical leakage follows the predicted dependence on regularisation, posterior variance and the number of released posterior sample-paths. Utility experiments on downstream posterior-sampling tasks identify noisy-observation regimes where privacycompatible regularisation preserves useful decisions with modest utility loss. When stronger privacy is needed, the intrinsic guarantee can be sharpened by adding calibrated GP noise, providing an explicit additional privacy knob.

1

INTRODUCTION

Protecting individual training examples is a central concern in modern machine learning. Released models can leak membership, memorised training fragments, or partial reconstructions of sensitive examples and attributes (Liu et al., 2023; Carlini et al., 2022, 2021; Balle et al., 2022; Haim et al., 2022; Ateniese et al., 2015). Standard DP methods mitigate this by injecting randomness during training or release (Dwork et al., 2006; Mironov, 2017; Hall et al., 2013; Abadi et al., 2016). However, recent work shows that randomness and regularisation already present in learning

algorithms can also shape privacy leakage: regularisation can reduce membership-inference risk (Tan et al., 2023), initialisation affects privacy–utility tradeoffs in overparameterised networks (Ye et al., 2023), and reconstruction attacks depend strongly on initialisation and regularisation (Haim et al., 2022; Balle et al., 2022; Buzaglo et al., 2023). This suggests that intrinsic model randomness may itself inhibit training-data inference even without explicit external DP noise. We study this effect in Gaussian-process (GP) posterior sampling, a tractable setting where randomness is part of the model rather than added externally. The mechanism maps a training set D to a posterior draw fD ∼ ΠD = GP(µD , σ 2 kD ) with the standard GP posterior mean and covariance, see Rasmussen and Williams (2005). The draw fD is not meant to be a private surrogate for the deterministic posterior mean. Rather, the posterior sample path is the released Bayesian object itself as in workflows based on simulation, uncertainty quantification, Thompson sampling and posterior-functional estimation. This gives a canonical setting for analysing how posterior randomness and regularisation jointly limit extraction of training-data information. The random mechanism D 7→ fD is a natural Bayesian release rule whose privacy guarantees are particularly relevant for one-shot posterior release, for release of multiple posterior draws via composition, and for any downstream object computed from such draws by post-processing. In particular, any output of a Bayesian algorithm that uses the private data only through posterior draws fD (e.g. the selected action in Thompson sampling or publicly released pollution map in environmental modeling) inherits the privacy guarantees we establish in this work. Posterior function draws are standard computational objects in Bayesian optimisation and Thompson sampling (Russo and Van Roy, 2014; Hernández-Lobato et al., 2017) including parallelised variants (Kandasamy et al., 2018). Posterior draws are also used in Bayesian workflows beyond optimisation where one repeatedly simulates from the posterior distribution and evaluates discrepancies or other posterior predictive checks (Gelman

Differential Privacy of Gaussian Process Posterior Sampling

et al., 1996; Gabry et al., 2019). Posterior sample draws of latent Gaussian field models used in spatial statistics and environmental modelling for estimating excursion sets, exceedance probabilities or contour uncertainty given some noisy latent field observations are likewise basic inferential objects (Bolin and Lindgren, 2015; French and Hoeting, 2016; Diggle and Ribeiro, 2007). More broadly, related privacy considerations may arise in distributed Bayesian computation where posterior draws are communicated across parties (Scott et al., 2016; Neiswanger et al., 2014). Threat Model and Main Contributions The adversary observes one or more released posterior sample paths fD drawn from the GP posterior conditioned on a private training set D. The release is function-valued i.e., the adversary may evaluate the released paths at arbitrary inputs. The posterior mean µD and covariance kD are not directly released. Our privacy guarantee does not rely on hiding the model assumptions – the adversary may also know the GP model, training data domain, and the size of D. 1. We derive the first explicit DP bounds for GP posterior sample-path release (Theorem 3.1). The bounds separate posterior-mean leakage from data-dependent covariance leakage, showing that posterior scale alone is insufficient and that effective ridge regularisation is essential for intrinsic privacy. 2. We test experimentally the theorem’s structural predictions. Our membership-inference attack shows that MIA success tracks the predicted dependence on r and σ. Utility experiments further show that in noisy-observation regimes privacy-compatible regularisation can preserve useful downstream decisions with modest utility loss, because the utility-favourable posterior is already substantially regularised. We state our bounds in terms of Rényi-DP guarantees for a single posterior draw which can be conveniently translated to bounds for multiple posterior draws and then translated to the approximate-DP bounds (i.e., (ε, δ)-DP) by the standard formulae, see Mironov (2017). Our DP bounds identify the kernel scale σ and the effective ridge r as the key privacy-controlling quantities and they show how posterior uncertainty and regularisation hinder extraction of training-data information. For posterior sample-path release, utility is inherently downstream-task dependent: the released function is useful only through post-processing. This is especially relevant in sensing applications where the measurement locations are themselves sensitive and the goal

is to release a spatial-based decision without revealing sensor locations. We study this through excursionset estimation in Section 5. The experiments show a clear privacy-utility transition: in noisy-observation regimes, the utility-optimal posterior is already substantially regularised, placing the mechanism closer to the strong-privacy region of our bounds. In these regimes meaningful intrinsic privacy comes with modest downstream utility loss, whereas low-noise regimes expose a sharper cost of private calibration. The tradeoff effect is less pronounced for larger training sets. Prior Work Posterior sampling has been connected to privacy in broader Bayesian settings. Wang et al. (2015) show that releasing one posterior sample is differentially private “for free” under bounded loglikelihood assumptions and Dimitrakakis et al. (2017) develop a general Bayesian posterior-sampling framework for DP. These results do not directly yield guarantees for the standard GP posterior-sampling mechanism considered here. The bounded log-likelihood condition of Wang et al. (2015, Theorem 1) fails for usual GP regression models because the latent field is unbounded. Similarly, Assumption 1 of Dimitrakakis et al. (2017) requires uniform Lipschitz continuity of the log-likelihood with respect to the data over all parameter values. In the GP setting these parameters are latent functions, so this condition is not satisfied. Geumlek et al. (2017) derive RDP guarantees for posterior sampling in exponential-family models, but GP posterior sampling does not have the conjugate exponential-family sufficient-statistic form required there. Prior GP-specific privacy work also differs in both mechanism and privacy model. Smith et al. (2018, 2021) assume public covariates X and private responses y, and privatise GP predictions by adding calibrated noise to posterior-mean predictions. In contrast, we treat the entire training set (covariates and responses) as private and analyse the release of a draw from the data-dependent GP posterior itself. To the best of our knowledge, these are the first explicit DP bounds for GP posterior sample-path release.

2

PROBLEM SETUP

We consider the setup where one aims to model an unknown function f∗ using a mean-zero GP prior with the covariance function σ 2 k : ΩX × ΩX → R where ΩX is the covariate domain and σ > 0. Given noisy observations y = (y1 , . . . , yn )T at covariates X = (x1 , . . . , xn ) we write D = (X, y). When conditioned on D, the posterior GP distribution is

Tomasz Maciazek

ΠD = GP(µD , σ 2 kD ) where  −1 (r) µD (x) = kX (x)T KXX y, r2 := λ2 /σ 2  −1 (r) kD (x, x′ ) = k (x, x′ ) − kX (x)T KXX kX (x′ ) (r)

with KXX := K(X, X) + r2 I, kX (x) := K(X, x). Here, λ2 is the noise variance assumed in the GP model, σ 2 is the kernel scale and r2 is the effective ridge sometimes referred to as the noise-to-signal variance ratio. Importantly, changing σ while keeping r fixed does not change the predictive mean µD – it only rescales the posterior covariance. We separate out the kernel scale σ 2 since it is the parameter that directly controls the posterior randomness. When k has bounded diagonal we rescale k so that supx∈ΩX k(x, x) = 1 in which case σ 2 has the standard interpretation as the prior variance scale. We denote by D ∼ D′ the replace-one neighbouring relation between datasets. Consider the random mechanism that assigns to D a single sample path from the GP -posterior  M : D 7→ fD ∼ GP µD , σ 2 kD . (1) In this paper we determine approximate-DP guarantees of this mechanism. Since M(D) is functionvalued, we use the formalism of differential privacy for function-valued mechanisms developed in Hall et al. (2013). More precisely, we view fD as a random el ement of the measurable space (F, A) = RΩX , C , where RΩX is the space of all real-valued functions on ΩX , and C is the canonical cylinder σ-algebra which makes all evaluation maps f 7→ f (x) measurable. The mechanism M is (ε, δ)-DP if for all D ∼ D′ and all S∈C Pr [fD ∈ S] ≤ eε Pr [fD′ ∈ S] + δ. (2) Crucially, it is enough to verify inequality (2) on finite evaluations, namely that Pr [fD (XT ) ∈ S] ≤ eε Pr [fD′ (XT ) ∈ S] + δ for every finite test set XT and every S being a measurable subset of R|XT | (see Hall et al., 2013, Proposition 6). To prove the (ϵ, δ)-DP guarantee for finite sample-path evaluations we turn to Rényi-DP (RDP) machinery. The mechanism D 7→ fD (XT ) is (α, εα )RDP (Mironov, 2017) iff for all D ∼ D′ we have Dα (ΠD (XT )∥ΠD′ (XT )) ≤ εα where ΠD (XT ) is the GP posterior distribution of fD (XT ) and Dα is the Rényi divergence of order α > 1. Having proved (α, εα )-RDP uniformly over all finite test sets XT , we apply the standard RDP-to-DP conversion which says that the mechanism is (ε′α , δ)-DP

1/δ with ε′α = εα + log α−1 (see Mironov, 2017, Proposition 3) for any 0 < δ < 1.

The RDP guarantees are also convenient when considering the mechanism that releases  L independently  (1) (L) sampled posterior draws ML : D 7→ fD , . . . , fD . If a single posterior draw satisfies (α, εα )-RDP uniformly over finite test sets, then by additive composition of RDP (Mironov, 2017, Proposition 1) ML satisfies (α, Lεα )-RDP. Consequently, ML is   log(1/δ) Lεα + , δ -DP (3) α−1 for full-function release. This accounting applies mutatis mutandis to adaptive releases where the round-ℓ posterior-sampling training set may be chosen as a function of previous draws and their postprocessed outputs, as in Thompson sampling. The final approximate-DP guarantee is obtained by optimising the converted bound over the admissible orders α > 1 for which the bound is finite. This optimisation can lead to sub-linear growth of the final approximateDP bound in L, see Section 4. The full function-valued release guarantee also permits lazy evaluation where the mechanism can answer fD (X1 ) and later sample fD (X2 ) | fD (X1 ) for an adaptively chosen X2 . The output is the same as evaluation of one fixed posterior draw, so no extra privacy budget is consumed.

3

RÉNYI-DP BOUNDS FOR GP POSTERIOR SAMPLING

As in standard additive-noise DP mechanisms the privacy guarantees depend on a problem-specific sensitivity of the non-private output which must be bounded for the release rule at hand (see e.g. Hall et al., 2013, Proposition 3). However, for the GP posteriorsampling mechanism, the released object is not a deterministic posterior mean plus data-independent external noise, but a draw from the full data-dependent posterior distribution. As we show, the relevant quantities are then the posterior variance scale and the common-core posterior mean sensitivity defined below. Definition 3.1. Let D ∼ D′ and let D− = D ∩ D′ denote the common core and let X− be the corresponding covariate set. Let kD− be the posterior covariance kernel obtained after conditioning on D− , as defined in Section 2. We define the posterior variance scale by Vn (r) :=

sup

sup kD− (x, x),

|X− |=n−1 x∈ΩX

and the common-core posterior mean sensitivity by ∆n (r) := sup ∥µD − µD′ ∥H− , D∼D ′

Differential Privacy of Gaussian Process Posterior Sampling

where H− denotes the reproducing kernel Hilbert space (RKHS) associated with the kernel kD− . For the RDP bounds presented in Theorem 3.1 below it is crucial to upper-bound Vn (r) and ∆n (r) as tightly as possible. Bounding Vn (r). In many cases Vn (r) can be conveniently upper-bounded or even computed exactly in closed form. When k(x, x) ≤ 1 for all x, then Vn (r) ≤ 1. If additionally k(x, x′ ) ≥ 0 for all x, x′ then n−1 , Vn (r) ≤ V̄n (r) := 1 − κ2 n − 1 + r2 where κ = inf x,x′ ∈ΩX k(x, x′ ). If on top of that k(x, x) = 1 for all x, then Vn (r) = V̄n (r), see Proposition 3 in Appendix B. This applies, for example, to standard normalised Matérn and RBF kernels. Bounding ∆n (r). The problem of finding bounds for ∆n (r) is more subtle. Let us restrict our considerations to kernels with bounded diagonal and responses bounded in absolute value by a constant MY > 0. Then, we have the generic bound (see Lemma B.1) √   √ MY n−1 ∆n (r) ≤ 1+ = O( n). r r Without additional information about the kernel k one √ cannot improve the O( n)-growth of ∆n (r). Indeed, in Lemma B.2 we give an example of a kernel (constructed via a particularly designed feature map) for √ which ∆n grows exactly at the n-rate. On the other hand, there exist typical situations where ∆n (r) is upper-bounded by a constant or even decays with n. For example, in Lemma B.4 we show that for the √ constant kernel k(x, x) ≡ 1 one has ∆n (r) = O(1/ n). More generally, if the kernel k is strictly positivevalued and has constant diagonal, then the existence of two distinct admissible datapoints with responses of opposite signs rules out any decay to zero of ∆n (r) with n. Indeed, Lemma B.3 shows that under these conditions ∆n (r) is bounded away from zero. Thus, the best one can hope for in this setting is ∆n (r) = O(1). In Lemmas B.7 and B.5 we show that the matching upper bounds ∆n (r) = O(1) hold for the one-dimensional exponential kernel and for the purely diagonal kernel in arbitrary √ dimension where ∆n (r) ≤ 2MY /r and ∆n (r) ≤ 2 2MY /r respectively. A canonical setting where a O(1)-bound holds for bounded-diagonal kernels is when the responses are generated by a fixed RKHS function, namely yi = f∗ (xi ) with f∗ ∈ Hk , see Lemma B.8. Then, ∆n (r) ≤ 2∥f∗ ∥Hk Vn (r)/(r2 + Vn (r)). The above examples illustrate how the abstract sensitivity quantities can be controlled in concrete regimes.

They are not additional assumptions of Theorem 3.1 which applies completely generally once valid bounds on Vn (r) and ∆n (r) are available. Main RDP Bound. We are now ready to state the main RDP bound which is uniform over all finite test sets XT . Theorem 3.1. Let D, D′ be replace-one neighbouring datasets of size n and let ΠD (XT ) and ΠD′ (XT ) denote the distributions of the values of the released posterior sample path on a test set XT . Then for every finite XT and every 1 < α < 1 + r2 /Vn (r) one has   Vn (r) Dα (ΠD (XT )∥ΠD′ (XT )) ≤ 2 ψα r2 (4)   2 α Vn (r) + r2 ∆n (r) + 2 r2 − (α − 1)Vn (r) σ where ψα (τ ) is a non-decreasing function of τ defined as (   1 + ατ 1 1 ψα (τ ) = max log(1 + τ ) − log , 2 2(α − 1) 1+τ ) 1 1 log (1 − τ (α − 1)) . − log(1 + τ ) − 2 2(α − 1) cn (r) In practice, we work with upper bounds Vn (r) ≤ V cn (r). Substituting these quantities for and ∆n (r) ≤ ∆ Vn (r) and ∆n (r) in (4) yields a valid RDP bound on cn (r). the admissible range 1 < α < 1 + r2 /V We prove Theorem 3.1 in Appendix A. The bound (4) separates the roles of the posterior scale σ and the effective ridge r. Increasing σ only attenuates the mean-sensitivity term through ∆n (r)/σ. It does not affect the covariance term 2ψα (Vn (r)/r2 ) which captures leakage from the data-dependent posterior covariance. This distinguishes posterior sampling from standard DP mechanisms based on data-independent noise: even in the limit σ → ∞, the centred posterior draw can leak through kD . Thus, meaningful intrinsic privacy also requires sufficient regularisation. Increasing r reduces Vn (r)/r2 , enlarges the admissible range of α, and decreases ∆n (r). This behaviour is visible in Figure 1a–b and is confirmed empirically in Section 4 where the membership-inference attack remains nontrivial at large σ but collapses toward random guessing as r increases. In Appendix A we also show that one can further enhance the intrinsic DP guarantees by adding noise to the released sample path via the mechanism D 7→ fD + g,

g ∼ GP(0, η 2 k)

(5)

Tomasz Maciazek

Figure 1: Examples of (ϵ, δ)-DP bounds constructed from Theorem 3.1 (log-log plots). We use Vn (r) = V̄n (r), MY = 1, n = 103 , δ = 5 · 10−4 , κ2 = 0.1. Left: the generic DP bound for ∆n (r) from Lemma B.1. Middle: more specialised setups lead to tighter DP bounds. It uses the tighter bounds in the RKHS response model from Lemma B.8 with ∥f∗ ∥Hk = 3. Right: the effect of using additional GP noise according to (5) with η = σ. This reduces ε by a factor 1 − 10 and the effect is strongest in less regularised regimes. where g is an independent sample path from the GP prior, see Appendix A.1. By increasing η one can obtain stronger privacy guarantees at fixed σ, r, albeit typically at the cost of reduced utility. The subsequent sections present experimental results whose implementation can be found on https:// github.com/tmaciazek/gaussian_process_dp.

4

MEMBERSHIP-INFERENCE ATTACK VS. OUR DP-BOUNDS

We use a LiRA-style MIA (Carlini et al., 2022) as an operational probe of the privacy effects captured by our bounds. We show that attack success tracks the predicted dependence on σ and r where increasing σ alone leaves residual membership signal, whereas increasing the effective ridge drives the attack performance to random guessing. The experiment is intended as a test of whether empirical leakage follows the dependence predicted by our bounds. This is not a comprehensive evaluation of MI risk. We assume the attacker can sample data from its prior distribution. We consider the one-dimensional setup where the responses are generated as yi = fstep (xi ), ΩX = [0, 1] and we sample the covariates from the uniform distribution. The step-function is fstep (x) = −1 if x < 1/2 and fstep (x) = 1 otherwise. We use the exponential kernel k(x, x′ ) = exp(−|x − x′ |/ℓ) for GP regression. The task is to decide whether the training set D (with |D| = 10) contained the point z0 := (x0 , y0 ) = (1/2, 1) given L released posterior (1) (L) sample paths fD , . . . , fD . To this end, we apply a

version of LiRA which uses the observables 1 fˆD (x0 ) := L v̂D (x0 ) :=

1 Lσ 2

L  X

L X

(l)

fD (x0 ),

l=1

2 (l) fD (x0 ) − fˆD (x0 ) .

l=1

Namely, we construct in- and out-distributions of   ˆ fD (x0 ), v̂(x0 ) by sampling training sets D that contain z0 and D that do not contain z0 . We subsequently empirically estimate the density of the in- and outdistributions (denoted by ρin and ρout respectively) and then apply the Neyman-Pearson binary hypothesis testing criterion   ρin fˆD (x0 ), v̂D (x0 )   > C, C ∈ R log ρout fˆD (x0 ), v̂D (x0 ) to decide whether the observed fˆD , v̂D came from D containing z0 . See Carlini et al. (2022) for more details of the LiRA attack. We use a latent Gaussian mixture model to estimate ρin/out , see Appendix C. Our choice of the noiseless regression function and the point z0 is deliberately favourable to the adversary as probing the data around the jump location strongly affects the GP posterior sample paths. Thus, including z0 creates a strong local membership signal. This is an appropriate stress-test for our DP claim which must control worst-case privacy leakage. We work with the 1D exponential kernel where by Proposition 3 we have Vn (r) = V̄n (r) with κ = exp(−1/ℓ) and Lemma B.7 gives us a specialised tighter sensitivity bound. For a fixed C in the Neyman-Pearson criterion we evaluate the true-positive rate (TPR) which is the fraction

Differential Privacy of Gaussian Process Posterior Sampling

Figure 2: MIA success and (ε, δ)-DP bounds for one released posterior sample path (L = 1, n = 10, δ = 0.05). Left: Excess TPR@FPR = 0.1 vs. ε for different σ. Error bands show standard deviation over 10 random seeds. Right: Privacy bound versus regularisation r. Dashed guides mark the ε = 10 threshold and the corresponding r values. Small ε is reached only with sufficient regularisation, even in the covarianceonly limit σ = ∞. of members correctly classified as members and the false-positive rate (FPR) which is the fraction of nonmembers incorrectly classified as members. The ROC curves are obtained by varying C. We report excess TPR@FPR = 10% i.e., T P R − 10% at F P R = 10%. Results at FPR = 1% are qualitatively similar and are deferred to Appendix C. Attack success vs. ε. Figure 2 assumes L = 1 and shows that the attack curves mirror our privacy bound. In the regime where the bound reaches single-digit ε, the LiRA-style attack is consistently driven close to random-guess performance, suggesting practical relevance of the bounds. This is in line with the recent observations (Lowy et al., 2024) that single-digit ε can suffice to defend against MIA under a realistic attacker model (such as ours) in which the adversary lacks exact knowledge of nearly all of the private data. When our bound remains large, large σ alone leaves residual membership signal, confirming that regularisation, not posterior scale alone, is crucial for good privacy. This is particularly evident when considering the formal limit σ → ∞ (with r-fixed) that corresponds to releasing covariance-only information via the centred posterior draw fD ∼ GP(0, kD ) where the RDP bound  simplifies to εα ≤ 2 ψα Vn (r)/r2 . Privacy composition with L. Figure 3 illustrates how MIA success scales with the number of released posterior sample paths L. As L becomes large, the attack approaches the non-private benchmark corresponding to releasing the exact posterior mean and variance. At the same time, the figure shows that optimising the RDP-to-DP conversion after composition yields a sub-linear scaling of the final (ε, δ)-DP bound with L despite the underlying RDP terms composing

Figure 3: Effect of releasing L posterior sample paths. The LiRA attack starts with random-guess effectiveness (T P R ≈ F P R) at L = 1 and strengthens with L approaching the non-private benchmark while the (ε, δ)-DP bound grows sub-linearly with L. Parameters: r = 1, σ = 5, n = 10, δ = 0.05. Error bands show standard deviation over 10 random seeds. additively as in Equation (3).

5

PRIVACY VS. UTILITY

We use synthetic experiments to isolate a privacyutility mechanism suggested by the privacy bounds. The goal is to show that private posterior sample-path releases can remain useful in regimes where our privacy bound is stronger. We use a simple one-dimensional model for which the privacy accounting is tight. This setting isolates the following transition: as observation noise makes stronger regularisation natural, posterior sampling moves into a stronger-privacy region while preserving high downstream utility. The experiment is intended to illustrate the possible usefulness of the studied DP mechanism. Assessing performance across broader tasks and more complex modelling settings requires further study. We consider a standard GP-regression setting where an unknown latent random field f∗ is observed only through noisy data D = (X, y), and inference is performed through a posterior GP, ΠD = GP(µD , σ 2 kD ) with hyperparameters θ = (σ, r, ℓ) corresponding to kernel scale, effective ridge and lengthscale. Abstractly, one would like to choose hyperparameters that minimise the average task loss, θ∗ ∈ arg min Ef∗ ,D [L(ΠD , f∗ )] . θ

(6)

Given θ one then evaluates the utility of the fitted posterior distribution U (ΠD , f∗ ). For ordinary regression the utility might be (negative) mean-squared error and the loss might be the predictive negative log-likelihood. Here we instead study a task for which posterior draws themselves are natural outputs: estimating excursion sets. Excursion-set problems are common in geospatial and environmental modelling Bolin and Lindgren (2015). For example, one may wish to publish a pollu-

Tomasz Maciazek

tion map from measurements collected at private sensor locations X. Given a threshold t ∈ R, the true excursion set is Ωt (f∗ ) := {x ∈ ΩX : f∗ (x) ≥ t}. In our experiments we take t = 0. The posterior assigns to each location x the excursion probability ! µD (x) − t p pD (x) := Pf ∼ΠD [f (x) ≥ t] = Φ , σ kD (x, x) where Φ is the CDF of normal distribution. Since the target is binary, we use the integrated binary crossentropy ℓBCE (s, p) := −s log p − (1 − s) log(1 − p) as the loss function: Z LBCE (ΠD , f∗ ) = ℓBCE (s∗ (x), pD (x)) dx ΩX

where s∗ (x) = 1 when f∗ (x) ≥ t and s∗ (x) = 0 otherwise. In the experiments this integral is approximated on a fixed mesh. Our released object is the excursion b The utility is the set-recovery accuracy set estimate Ω. which we measure by intersection-over-union (IoU),

where |Ω| denotes the volume of Ω. The non-private benchmark for hyperparameters θ∗ is b D (θ∗ ) = {x ∈ ΩX : pD (x) ≥ C}, Ω where the cutoff C is chosen by validation. The private release uses hyperparameters selected by the constrained search min θ: εL (θ)≤ε0

Ef∗ ,D [L(ΠD , f∗ )] ,

0.1

0.3

0.5

0.6

NSR 0.11 0.43 1.00 1.50 ε, unconstrained 2· 104 500 63.6 21.4 deff , unconstrained 52 17.7 8.47 6.01 deff , ε < 10 2.26 2.21 2.15 2.08 relative BCE increase 115% 33.5% 8.8% 2.2% IoU, unconstrained 0.957 0.922 0.874 0.844 1-path IoU, ε < 10 0.854 0.833 0.801 0.772 relative IoU gap 10.1% 8.7% 7.7% 7.3%

Table 1: Median over 103 draws of (D, f∗ )-pairs across a range NSR values. The private candidate is the best average-BCE-selected setting satisfying ε < 10. The 1-path IoU is averaged over 50 posterior draws. (L)

b . Since this is a postthe GP posterior and release Ω D processing of the sampled paths, this set-valued release b (L) inherits the same (εL , δ) privacy guarantee as the Ω D L released posterior paths. The connection to randomised decision-making is conceptual rather than a direct application of the results of Dimitrakakis et al. (2017) – their assumptions do not cover the usual GP regression posterior over unbounded latent functions which is why we derive GP-specific bounds. Empirical results. We generate the responses by sampling f∗ ∼ GP(0, k̃) given noisy observations generated according to the following bounded-response model

b ∩ Ωt (f∗ )  Ω b , Uθ (ΠD , f∗ ) = IoU Ω, Ωt (f∗ ) := b ∪ Ωt (f∗ ) Ω 

priv θL ∈ arg

(7)

where εL (θ) is the (ε, δ)-DP privacy bound for releasing L independent posterior sample paths and ε0 is the target privacy budget. Given these hyperparameters, the released randomized estimate is ( ) L n o X 1 (L) (ℓ) b Ω x ∈ ΩX : 1 fD (x) ≥ t ≥ c , D = L ℓ=1

(ℓ)

where c is fixed by validation and fD are i.i.d. from priv the GP posterior using θ = θL . This operational approach is in the spirit of the randomized decision-theoretic view of posterior sampling developed by Dimitrakakis et al. (2017). We use posterior draws as a randomized answer to a down(1) (L) stream excursion-set query: draw fD , . . . , fD from

yi = (1 − Mξ )

f∗ (xi ) + ξi , ∥f∗ ∥∞

i.i.d.

ξi ∼ Pξ

(8)

with Pξ = Unif[−Mξ , Mξ ], 0 < Mξ < 1 and ∥f∗ ∥∞ := supx∈ΩX |f∗ (x)|. In the experiment we take k̃(x, x′ ) = exp(−|x − x′ |) on ΩX = [0, 1] and use the exponential GP-kernel for regression. Such a GP model allows us to apply Lemma B.7 to find tighter DP-bounds. The GP-regression hyperparameters are θ = (r, σ, ℓ). Such response model has amplitude noise-to-signal ratio N SR := Mξ /(1 − Mξ ). We work with L = 1, |D| = 100 and the privacy budget ε0 = 10 with priv δ = 0.005. The parameters θ∗ and θL are selected via several rounds of grid-search minimisation according to Equations (6) and (7) on a fixed sample of 103 pairs (D, f∗ ). The thresholds c, C are chosen via another grid-search on a separate validation set by optimising the expected IoU. Table 1 shows the transition across noise levels. As Mξ increases the unconstrained utility optimum becomes naturally more regularised: both ε and deff decrease sharply. Thus, the cost of imposing ε < 10 also decreases with the relative BCE increase falling from 115% at NSR = 0.11 to 2.2% at NSR = 1.50. The intermediate case NSR = 1 illustrates the typical tradeoff where the privacy constraint reduces deff from 8.47

Differential Privacy of Gaussian Process Posterior Sampling

6

DISCUSSION AND LIMITATIONS

Our DP bounds assume a fixed kernel and fixed hyperparameters, including r, σ and ℓ. They apply directly when these choices are independent of the private dataset, e.g. selected using public or non-private set-aside data. If hyperparameters are selected using the private training set and then used or reported, that selection must be privacy-accounted. Any RDP guarantee for this step can be combined with ours by RDP composition, e.g. for hyperparameters drawn from a privacy-accounted posterior in hierarchical GP.

Figure 4: Greater London 2018 house-price excursion experiment. The task is to estimate the region where the median log sale price exceeds 13 (about £442k). The private posterior draws (dashed lines) retain the dominant high-price region in central London, but typically miss or smooth out smaller high-price islands outside central London illustrating the additional randomisation and regularisation induced by the privacy-constrained posterior. ε = 10, δ = n−1.1 with n = 1075 hexagons.

to 2.15 while increasing median BCE by only 8.8%. The released object is the thresholded posterior sample path itself and by post-processing it inherits the onepath privacy guarantee. Its median IoU remains high, e.g. 0.8 at NSR = 1, only 7.7% below the non-private posterior-probability benchmark. The IoU gap is more stable than the BCE gap because IoU measures thresholded geometry rather than posterior calibration. In low-noise regimes privacy is more costly because the unconstrained optimum fits the data closely and lies far from the privacy-feasible region. At the same time, the excursion problem is easier because the observations are more informative. These effects partially offset each other: the soft posterior field pD (x) may degrade while the released excursion set remains accurate. For more details see Appendix D. Table 3 in Appendix D shows that the utility tradeoff becomes smaller when |D| grows at fixed NSR. Appendix D.1 presents a 2D excursion set study with excursion set visualisations. As a real-data illustration, we apply the same procedure to HM Land Registry (2026) for 2018 Greater London leasehold-flat transactions. Transactions are aggregated onto a hexagonal lattice and the GP mechanism is applied to the clipped and rescaled centred responses, see Appendix D.2 for details. Figure 4 shows the resulting non-private posterior excursion boundary and private one-path excursion boundaries.

Meaningful small-ε guarantees require moderate-tostrong effective regularisation which may be restrictive for membership privacy. Larger ε values may be relevant for higher-impact threats such as full-example reconstruction: prior work on private neural networks shows that reconstruction can be mitigated even at ε ∼ 102 − 103 (Balle et al., 2022; Hayes et al., 2023). Whether similar phenomenon holds for GP posterior sampling is open and requires developing new tailored reconstruction attacks. Our bounds also require ∆n (r) < ∞ which relies on bounded responses. This is an issue which is ubiquitous in DP analyses – unbounded responses can lead to unbounded sensitivity. Standard mitigations apply a training-data preprocessing step such as response clipping while accounting for the resulting utility tradeoff. Theorem 3.1 is general, but the practical strength of the bound depends on how tightly one can control ∆n (r). The generic bound (Lemma B.1) can be conservative and deriving sharper ∆n (r) bounds for other kernels and response models beyond those developed here is an important direction for future work. Conclusions We show that GP posterior sampling can provide intrinsic privacy guarantees in regimes controlled by the posterior covariance and effective regularisation. Our RDP bounds make this explicit and extend to multiple posterior draws by composition. Empirically, membership-inference attacks follow the predicted dependence on r, σ, and the number of released paths L. Interestingly, covariance-only privacy leakage is possible even when centered posterior sample paths fD ∼ GP(0, kD ) are released. The utility experiments show that for randomized excursionset estimation privacy-compatible regularisation can preserve utility even when it degrades posterior calibration. In particular, noisy regimes in which regularisation is natural yield meaningful single-path privacy guarantees with useful released excursion sets. Thus posterior randomness can be a genuine privacy resource, but not an unconditional one – useful intrinsic privacy requires enough regularisation to control both posterior-mean and posterior-covariance leakage.

Tomasz Maciazek

Acknowledgements TM has been fully funded by His Majesty’s Government. This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol - http://www.bristol. ac.uk/acrc/. References Martin Abadi, Andy Chu, Ian Goodfellow, H. Brendan McMahan, Ilya Mironov, Kunal Talwar, and Li Zhang. Deep learning with differential privacy. In Proceedings of the 2016 ACM SIGSAC Conference on Computer and Communications Security, 2016. Giuseppe Ateniese, Luigi V. Mancini, Angelo Spognardi, Antonio Villani, Domenico Vitali, and Giovanni Felici. Hacking smart machines with smarter ones: How to extract meaningful data from machine learning classifiers. Int. J. Secur. Netw., 10 (3):137–150, September 2015. ISSN 1747-8405. doi: 10.1504/IJSN.2015.071829. Borja Balle, Giovanni Cherubin, and Jamie Hayes. Reconstructing training data with informed adversaries. In 2022 IEEE Symposium on Security and Privacy (SP), pages 1138–1156, 2022. doi: 10.1109/ SP46214.2022.9833677. David Bolin and Finn Lindgren. Excursion and contour uncertainty regions for latent gaussian models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 77(1):85–106, 2015. doi: https://doi.org/10.1111/rssb.12055. Gon Buzaglo, Niv Haim, Gilad Yehudai, Gal Vardi, Yakir Oz, Yaniv Nikankin, and Michal Irani. Deconstructing data reconstruction: Multiclass, weight decay and general losses. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 51515–51535. Curran Associates, Inc., 2023. Nicholas Carlini, Florian Tramèr, Eric Wallace, Matthew Jagielski, Ariel Herbert-Voss, Katherine Lee, Adam Roberts, Tom Brown, Dawn Song, Úlfar Erlingsson, Alina Oprea, and Colin Raffel. Extracting training data from large language models. In 30th USENIX Security Symposium (USENIX Security 21), pages 2633–2650. USENIX Association, August 2021. ISBN 978-1-939133-24-3. Nicholas Carlini, Steve Chien, Milad Nasr, Shuang Song, Andreas Terzis, and Florian Tramèr. Membership inference attacks from first principles. In 2022 IEEE Symposium on Security and Privacy (SP), pages 1897–1914, 2022. doi: 10.1109/SP46214.2022. 9833649.

Peter J. Diggle and Paulo J. Ribeiro. Model-based Geostatistics. Springer Series in Statistics. Springer, New York, NY, 1 edition, 2007. ISBN 978-0-38732907-9. doi: 10.1007/978-0-387-48536-2. Christos Dimitrakakis, Blaine Nelson, Zuhe Zhang, Aikaterini Mitrokotsa, and Benjamin I. P. Rubinstein. Differential privacy for bayesian inference through posterior sampling. Journal of Machine Learning Research, 18(11):1–39, 2017. Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam Smith. Calibrating noise to sensitivity in private data analysis. In Shai Halevi and Tal Rabin, editors, Theory of Cryptography, pages 265– 284, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg. ISBN 978-3-540-32732-5. Joshua P. French and Jennifer A. Hoeting. Credible regions for exceedance sets of geostatistical data. Environmetrics, 27(1):4–14, 2016. doi: https://doi. org/10.1002/env.2371. Jonah Gabry, Daniel Simpson, Aki Vehtari, Michael Betancourt, and Andrew Gelman. Visualization in bayesian workflow. Journal of the Royal Statistical Society Series A: Statistics in Society, 182(2):389– 402, 02 2019. ISSN 0964-1998. doi: 10.1111/rssa. 12378. Andrew Gelman, Xiao-Li Meng, and Hal Stern. Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica, 6(4):733–760, 1996. ISSN 10170405, 19968507. Joseph Geumlek, Shuang Song, and Kamalika Chaudhuri. Rényi differential privacy mechanisms for posterior sampling. In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, page 5295–5304, Red Hook, NY, USA, 2017. Curran Associates Inc. ISBN 9781510860964. M. Gil, F. Alajaji, and T. Linder. Rényi divergence measures for commonly used univariate continuous distributions. Information Sciences, 249:124–131, 2013. ISSN 0020-0255. doi: https://doi.org/10. 1016/j.ins.2013.06.018. Niv Haim, Gal Vardi, Gilad Yehudai, Ohad Shamir, and Michal Irani. Reconstructing training data from trained neural networks. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 22911–22924. Curran Associates, Inc., 2022. Rob Hall, Alessandro Rinaldo, and Larry Wasserman. Differential privacy for functions and functional data. J. Mach. Learn. Res., 14(1):703–727, February 2013. ISSN 1532-4435.

Differential Privacy of Gaussian Process Posterior Sampling

Jamie Hayes, Borja Balle, and Saeed Mahloujifar. Bounding training data reconstruction in DP-SGD. In Thirty-seventh Conference on Neural Information Processing Systems, 2023. URL https:// openreview.net/forum?id=7LZ4tZrYlx. José Miguel Hernández-Lobato, James Requeima, Edward O. Pyzer-Knapp, and Alán Aspuru-Guzik. Parallel and distributed thompson sampling for large-scale accelerated exploration of chemical space. In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, page 1470–1479. JMLR.org, 2017. HM Land Registry. Price Paid Data. https://www. gov.uk/government/statistical-data-sets/ price-paid-data-downloads, 2026. Contains HM Land Registry data © Crown copyright and database right 2021. Licensed under the Open Government Licence v3.0. Accessed 2026-06-10. Kirthevasan Kandasamy, Akshay Krishnamurthy, Jeff Schneider, and Barnabas Poczos. Parallelised bayesian optimisation via thompson sampling. In Amos Storkey and Fernando Perez-Cruz, editors, Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pages 133–142. PMLR, 09–11 Apr 2018.

Daniel Russo and Benjamin Van Roy. Learning to optimize via posterior sampling. Math. Oper. Res., 39(4):1221–1243, November 2014. ISSN 0364-765X. doi: 10.1287/moor.2014.0650. Steven L. Scott, Alexander W. Blocker, Fernando V. Bonassi, Hugh A. Chipman, Edward I. George, and Robert E. McCulloch. Bayes and Big Data: The Consensus Monte Carlo Algorithm. International Journal of Management Science and Engineering Management, 11:78–88, 2016. Michael T. Smith, Mauricio A. Álvarez, Max Zwiessele, and Neil D. Lawrence. Differentially private regression with Gaussian processes. In Amos Storkey and Fernando Perez-Cruz, editors, Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pages 1195– 1203. PMLR, 09–11 Apr 2018. Michael Thomas Smith, Mauricio A. Álvarez, and Neil D. Lawrence. Differentially private regression and classification with sparse gaussian processes. J. Mach. Learn. Res., 22(1), January 2021. ISSN 15324435.

Lan Liu, Yi Wang, Gaoyang Liu, Kai Peng, and Chen Wang. Membership inference attacks against machine learning models via prediction sensitivity. IEEE Transactions on Dependable and Secure Computing, 20(3):2341–2347, 2023. doi: 10.1109/TDSC. 2022.3180828.

Jasper Tan, Daniel LeJeune, Blake Mason, Hamid Javadi, and Richard G. Baraniuk. A blessing of dimensionality in membership inference through regularization. In Francisco Ruiz, Jennifer Dy, and Jan-Willem van de Meent, editors, Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, volume 206 of Proceedings of Machine Learning Research, pages 10968–10993. PMLR, 25–27 Apr 2023.

Andrew Lowy, Zhuohang Li, Jing Liu, Toshiaki KoikeAkino, Kieran Parsons, and Ye Wang. Why does differential privacy with large ε defend against practical membership inference attacks? In Proceedings of the 5th AAAI Workshop on Privacy-Preserving Artificial Intelligence (PPAI), 2024. URL https: //aaai-ppai24.github.io/#accepted_papers.

Yu-Xiang Wang, Stephen Fienberg, and Alex Smola. Privacy for free: Posterior sampling and stochastic gradient monte carlo. In Francis Bach and David Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 2493– 2502, Lille, France, 07–09 Jul 2015. PMLR.

Ilya Mironov. Rényi differential privacy. In 2017 IEEE 30th Computer Security Foundations Symposium (CSF), pages 263–275, 2017. doi: 10.1109/ CSF.2017.11.

Jiayuan Ye, Zhenyu Zhu, Fanghui Liu, Reza Shokri, and Volkan Cevher. Initialization matters: Privacyutility analysis of overparameterized neural networks. In Thirty-seventh Conference on Neural Information Processing Systems, 2023.

Willie Neiswanger, Chong Wang, and Eric P. Xing. Asymptotically exact, embarrassingly parallel mcmc. In Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, UAI’14, page 623–632, Arlington, Virginia, USA, 2014. AUAI Press. ISBN 9780974903910. Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 11 2005. ISBN 9780262256834. doi: 10.7551/mitpress/3206.001.0001.

Supplementary Materials

A

PROOF OF THE RDP BOUND

Throughout this section, let D = D− ∪ {(x, y)} and D′ = D− ∪ {(x′ , y ′ )} be replace-one neighbouring datasets, where D− := ((x1 , y1 ), . . . , (xn−1 , yn−1 )) is their common core. Write also X− := (x1 , . . . , xn−1 ) and y− := (y1 , . . . , yn−1 )T . Define the common-core posterior mean and covariance as −1 µD− (X) := K(X, X− ) K(X− , X− ) + r2 I y− , −1 K(X− , X ′ ). KD− (X, X ′ ) := K(X, X ′ ) − K(X, X− ) K(X− , X− ) + r2 I For a finite test set XT , set C := KD− (XT , XT ),

c(x) := kD− (XT , x),

v(x) := kD− (x, x).

The identities for the posterior mean and covariance µD (XT ) = µD− (XT ) +

y − µD− (x) c(x), v(x) + r2

y ′ − µD− (x′ ) c(x′ ), v(x′ ) + r2 c(x′ ).c(x′ )T KD′ = C − v(x′ ) + r2

µD′ (XT ) = µD− (XT ) + c(x).c(x)T , v(x) + r2

KD := kD (XT , XT ) = C −

(9)

are the standard one-step update formulas obtained by writing the full posterior expressions for D = D− ∪ {(x, y)} and D′ = D− ∪ {(x′ , y ′ )} from Section 2 in block form relative to the common core D− and taking the corresponding Schur complement and using the block-matrix inverse identity. We also define the common-core posterior variance scale VD− (r) and the uniform posterior variance scale Vn (r) as VD− (r) := sup kD− (x, x),

Vn (r) :=

x∈ΩX

sup

VD− (r).

|D− |=n−1

Since 0 ≤ kD− (x, x) ≤ k(x, x) and k(x, x) ≤ 1, we have 0 ≤ VD− (r), Vn (r) ≤ 1. Finally, for the posterior mean we define the uniform common-core sensitivity 1/2 ∆n (r) := sup sup δTT kD− (XT , XT )† δT , δT := µD (XT ) − µD′ (XT ), |XT |<∞ D∼D ′

where A† denotes the Moore–Penrose pseudoinverse. If kD− (XT , XT ) is invertible this reduces to the usual inverse. Lemma A.1 (Covariance sandwich). For every replace-one neighbouring pair D, D′ with common core D− and every finite test set XT r2 C ⪯ KD ⪯ C, VD− (r) + r2

r2 C ⪯ KD′ ⪯ C. VD− (r) + r2

Consequently, one has r2 KD′ ⪯ KD ⪯ VD− (r) + r2



VD− (r) 1+ r2

Moreover, for Kα := αKD′ + (1 − α)KD one has Kα ⪰

r2 − (α − 1)VD− (r) C. VD− (r) + r2 2

Thus, Kα is positive semidefinite whenever 1 < α < 1 + VDr (r) . −

 KD′ .

Differential Privacy of Gaussian Process Posterior Sampling

Proof. Consider first KD := kD (XT , XT ). Since  KD− (XT ∪ {x}, XT ∪ {x}) =

 c(x) v(x)

C c(x)T

is positive semidefinite, for every a, t we have     C c(x) a T a t = aT Ca + 2taT c(x) + t2 v(x) ≥ 0. t c(x)T v(x) The discriminant of this polynomial in t must be non-positive which is equivalent to aT c(x)

2

≤ v(x) aT Ca.

(10)

Equivalently, c(x).c(x)T ⪯ vC. Hence by Equations (9) KD = C −

v(x) r2 c(x).c(x)T ⪰ C − C = C, v(x) + r2 v(x) + r2 v(x) + r2

KD′ ⪰

r2 C, v(x′ ) + r2

where we applied the same argument to x′ to obtain the corresponding bound for KD′ . Since v(x′ ), v(x) ≤ VD− (r), this yields r2 r2 ′ ⪰ KD ⪰ C, K C, (11) D VD− (r) + r2 VD− (r) + r2 The upper bound KD ⪯ C is immediate from the identity for KD . Indeed, by Equations (9) C − KD =

c(x).c(x)T ⪰ 0. v(x) + r2

Combining this with (11) yields KD ⪯ C ⪯ (1 + τ− (r)) KD′ ,

τ− (r) :=

VD− (r) . r2

This proves KD ⪯ (1 + τ− (r))KD′ and interchanging D and D′ yields KD′ ⪯ (1 + τ− (r))KD which is equivalent to 1 1 KD′ ⪯ KD ⪯ (1 + τ− (r))KD′ , KD ⪯ KD′ ⪯ (1 + τ− (r))KD . 1 + τ− (r) 1 + τ− (r) Finally, using the lower bound for KD′ and the upper bound for KD , Kα = αKD′ − (α − 1)KD ⪰ α =

r2 C − (α − 1)C = VD− (r) + r2



 αr2 − (α − 1) C VD− (r) + r2

r2 − (α − 1)VD− (r) C. VD− (r) + r2

This is positive semidefinite whenever r2 − (α − 1)VD− (r) > 0. Proposition 1 (RDP bound). Let D, D′ be replace-one neighbouring labelled datasets of size n. Let ΠD (XT ) and ΠD′ (XT ) denote the (Gaussian) distributions of the released posterior sample vector on a finite test set XT . Then for every r2 1<α<1+ , Vn (r) one has  Dα (ΠD (XT )∥ΠD′ (XT )) ≤ 2 ψα where

 ψα (τ ) :=

max λ∈[ 1/(1+τ ), 1+τ ]

Vn (r) r2



 α Vn (r) + r2 + 2 2 ∆2 (r), 2σ (r − (α − 1)Vn (r)) n

 1 1 − log λ − log (α − (α − 1)λ) . 2 2(α − 1)

Tomasz Maciazek

Proof. Write ΠD (XT ) = N (µD , σ 2 KD ),

ΠD′ (XT ) = N (µD′ , σ 2 KD′ ).

For multivariate Gaussian measures the Rényi divergence is given by the following formula (see Gil et al., 2013, Table 2), Dα (ΠD (XT )∥ΠD′ (XT )) =

1 det(Kα ) α (µD − µD′ )T Kα−1 (µD − µD′ ) − , log 2σ 2 2(α − 1) det(KD )1−α det(KD′ )α

(12)

where Kα := αKD′ + (1 − α)KD ≻ 0. Let us first consider the finiteness condition in the Gaussian Rényi formula. The sandwich bounds from Lemma A.1 imply that Kα ≻ 0 whenever C ≻ 0 and 1 < α < 1 + r2 /Vn (r). In general, KD , KD′ , and Kα may be singular as matrices on R|XT | . However, the same sandwich bounds imply that on this range of α they have the same nullspace as C. From (9) µD (XT ) − µD′ (XT ) =

y ′ − µD− (x′ ) y − µD− (x) c(x) − c(x′ ). v(x) + r2 v(x′ ) + r2

We also have that c(x), c(x′ ) ∈ range(C). To see this, consider inequality (10) from Lemma A.1 which says that 2 2 for any vector a we have aT c(x) ≤ v(x) aT Ca. Taking a ∈ ker(C) gives aT c(x) ≤ 0, so we neccessarily have aT c(x) = 0 for any a ∈ ker(C). In other words, c(x) ∈ ker(C)⊥ for any x. Since C is symmetric we have ker(C)⊥ = range(C), thus c(x) ∈ range(C) for any x. Hence, µD (XT ) − µD′ (XT ) ∈ range(C). Thus, the Gaussian measures ΠD (XT ) and ΠD′ (XT ) are supported on the same affine subspace determined by C. So, the Gaussian Rényi formula applies and remains finite on this common support even when Kα is singular with ordinary inverses and determinants replaced by Moore-Penrose pseudoinverses and pseudodeterminants. Next, we treat the two terms in (12) separately. Let us start with the determinant term. By Lemma A.1, 1 KD′ ⪯ KD ⪯ (1 + τ− (r)) KD′ , 1 + τ− (r) −1/2

Next, conjugating the above inequality by KD′

V− (r) . r2

τ− (r) := −1/2

and noting that KD′

−1/2

KD′ KD′

= 1 yields

1 −1/2 1 ⪯ KD−1/2 KD KD′ ⪯ (1 + τ− (r)) 1. ′ 1 + τ− (r) −1/2

Define R := KD′

(13)

−1/2

KD KD′ 

. By (13), the eigenvalues of R all lie in the interval    1 1 Vn (r) , 1 + τ− (r) ⊂ , 1 + τn (r) , τn (r) := . 1 + τ− (r) 1 + τn (r) r2

We also have 1/2

1/2

Kα = KD′ (αI − (α − 1)R) KD′

(14)

Moreover, by (9) KD − K D ′ =

c(x′ ).c(x′ )T c(x).c(x)T − , ′ 2 v(x ) + r v(x) + r2

−1/2

−1/2

so rank(KD − KD′ ) ≤ 2. Note that R − I = KD′ (KD − KD′ )KD′ , thus we also have that rank(R − I) ≤ 2. Let λ1 , . . . , λm denote the eigenvalues of R. Using the identity (14), the determinant term becomes 1 1 det Kα det(αI − (α − 1)R) log =− log 1−α α 2(α − 1) (det KD ) (det KD′ ) 2(α − 1) (det R)1−α   m X 1 1 = − log λj − log (α − (α − 1)λj ) ≤ 2 ψα (τ− (r)) ≤ 2 ψα (τn (r)) , 2 2(α − 1) j=1

Differential Privacy of Gaussian Process Posterior Sampling

where in the penultimate inequality we have used the fact that rank(R − I) ≤ 2 which implies that at most 2 of the λj ’s differ from 1 and contribute to the sum. It remains to bound the mean-term in (12). By Lemma A.1, Kα ⪰

r2 − (α − 1)V− (r) C. V− (r) + r2

Since Kα and C have the same nullspace on the admissible range of α, we may invert this inequality on their common support by taking the Moore-Penrose pseudoinverse of both sides to obtain Kα† ⪯

V− (r) + r2 C †. r2 − (α − 1)V− (r)

Therefore, with δT := µD (XT ) − µD′ (XT ),  α V− (r) + r2 α T † δ K δT ≤ δ T C † δT . 2σ 2 T α 2σ 2 (r2 − (α − 1)V− (r)) T Taking the supremum over all neighbouring pairs and all finite test sets XT , and using V− (r) ≤ Vn (r), yields  α Vn (r) + r2 α T † δ K δT ≤ ∆2 (r). 2σ 2 T α 2σ 2 (r2 − (α − 1)Vn (r)) n Combining the determinant and mean bounds proves the claim. The following lemma allows us to explicitly take the supremum over XT in the definition of ∆n (r) as originally stated in Definition 3.1. Lemma A.2 (∆n (r) as RKHS norm). Let D = D− ∪{(x, y)} and D′ = D− ∪{(x′ , y ′ )} be replace-one neighbouring datasets. Then µD − µD′ ∈ H− , where H− denotes the reproducing kernel Hilbert space associated with the common-core posterior covariance kernel kD− . For any finite test set XT ⊂ ΩX define δT := µD (XT )−µD′ (XT ). Then δTT kD− (XT , XT )† δT ≤ ∥µD − µD′ ∥2H− . Equality holds whenever {x, x′ } ⊂ XT . Consequently, ∆n (r) = sup ∥µD − µD′ ∥H− . D∼D ′

Proof. Fix D, D′ and write g := µD − µD′ . By the one-step update formula (9), if x and x′ are the differing inputs in D and D′ , then g(·) = η kD− (·, x) − η ′ kD− (·, x′ ) for suitable scalars η, η ′ . Hence g ∈ H− . Next, we use the general standard fact that for any h from RKHS of some kernel k̃ and for any finite sequence of distinct points XT ⊂ ΩX we have (see Hall et al., 2013, Proposition 8) h(XT )T K̃(XT , XT )−1 h(XT ) ≤ ∥h∥2Hk̃ . (15) In the possibly singularPcase the same statement holds with the Moore-Penrose pseudoinverse. P The equality in m (15) holds when h(·) = i=1 ai k̃(·, xi ) with (x1 , . . . , xm ) ⊂ XT . To see this, expand ∥h∥2H = i,j ai aj k̃(xj , xi ). k̃

On the other hand, denoting G := K̃(XT , XT ) we have X X h(XT )T G̃† h(XT ) = ai aj k̃(XT , xj )T G† k̃(XT , xi ) = ai aj eTj GG† Gei i,j

=

X i,j

i,j

ai aj eTj Gei =

X

ai aj k̃(xj , xi ),

i,j

where in the second equality we have used the fact that k̃(XT , xi ) = Gei with ei being the i-th unit vector in R|XT | , in the third equality we have used the Moore-Penrose property GG† G = G and in the last equality we have used the fact that eTj Gei = Gji = k̃(xj , xi ).

Tomasz Maciazek

The final inequality is obtained by putting h = g and k̃ = kD− in (15). This means that sup δTT kD− (XT , XT )† δT = ∥µD − µD′ ∥2H− ,

|XT |<∞

proving that ∆n (r)2 = supD∼D′ ∥µD − µD′ ∥2H− . The following lemma proves the explicit form of ψα (τ ) stated in Theorem 3.1. Lemma A.3. Let τ ≥ 0 and 1 < α < 1 + 1/τ with the convention that 1/τ = +∞ when τ = 0. Define   1 1 ψα (τ ) := max − log λ − log (α − (α − 1)λ) . 2 2(α − 1) λ∈[ 1/(1+τ ), 1+τ ] Then,  ψα (τ ) = max

1 1 log(1 + τ ) − log 2 2(α − 1)



1 + ατ 1+τ



 1 1 , − log(1 + τ ) − log (1 − τ (α − 1)) . 2 2(α − 1)

Proof. It is straightforward to show that 1 1 log (α − (α − 1)λ) ϕα (λ) := − log λ − 2 2(α − 1) is convex on the interval [ 1/(1 + τ ), 1 + τ ] by calculating its second derivative. Thus, it’s maximum is attained at one of the interval’s endpoints. This yields the final result. A.1

Adding prior GP noise to the released sample path

We now consider a simple modification of the release mechanism in which one first draws  fD (XT ) ∼ N µD (XT ), σ 2 KD , KD := kD (XT , XT ), from the posterior and then adds an independent Gaussian perturbation drawn from the GP prior  gT ∼ N 0, η 2 K(XT , XT ) , η ≥ 0. The released vector is therefore feD (XT ) := fD (XT ) + gT . e D denotes the resulting sample path distribution, then If Π   e D (XT ) = N µD (XT ), K eD , e D := σ 2 KD + η 2 K(XT , XT ). Π K The next proposition shows that the same proof strategy as in Proposition 1 yields an RDP bound for this modified mechanism. The additional GP noise improves both the covariance term and the posterior-mean term of the original RDP bound. Proposition 2 (RDP bound with additional GP noise). Let D, D′ be replace-one neighbouring labelled datasets of size n. On a finite test set XT consider the modified release  feD (XT ) = fD (XT ) + gT , gT ∼ N 0, η 2 K(XT , XT ) , where gT is independent of fD (XT ). Denote the distribution of feD (XT ) by   e D (XT ) = N µD (XT ), K eD , e D := σ 2 KD + η 2 K(XT , XT ), Π K e D′ (XT ) analogously. Then, for every α such that and define Π 1<α<1+

1 , τen (r, η)

τen (r, η) :=

σ 2 Vn (r) 2 2 σ r + η 2 (Vn (r) + r2 )

one has   α Vn (r) + r2 e D (XT ) Π e D′ (XT ) ≤ 2 ψα (e Dα Π τn (r, η)) + ∆2 (r). 2 σ 2 (r2 − (α − 1)Vn (r)) + η 2 (Vn (r) + r2 ) n

Differential Privacy of Gaussian Process Posterior Sampling

Proof. Let KT := K(XT , XT ). Since C = KT − K(XT , X− ) K(X− , X− ) + r2 I

−1

K(X− , XT ),

we can write H ⪰ 0.

KT = C + H,

Using the one-step update covariance identities (9), we obtain ′

T

′ T

e D′ = η 2 H + (σ 2 + η 2 )C − σ 2 c(x ).c(x ) . K v(x′ ) + r2

e D = η 2 H + (σ 2 + η 2 )C − σ 2 c(x).c(x) , K v(x) + r2

(16)

As in the proof of Lemma A.1, the positive-semidefiniteness of block matrices involving C, c(x), v(x) and C, c(x), v(x′ ) implies that c(x).c(x)T ⪯ v(x) C and c(x′ ).c(x′ )T ⪯ v(x′ ) C. Since v(x), v(x′ ) ≤ V− (r), both noisy covariance matrices satisfy   r2 2 2 2 e e C. KD , KD′ ⪰ η H + η + σ V− (r) + r2 Since c(x).c(x)T , c(x′ ).c(x′ )T ⪰ 0, we also have eD, K e D′ ⪯ η 2 H + (σ 2 + η 2 )C. K Next, define τe− (r, η) through 1 + τe− (r, η) =

σ2 + η2 . 2 − (r) + r )

η 2 + σ 2 r2 /(V

The preceding bounds imply the sandwich relation 1 e D′ ⪯ K e D ⪯ (1 + τe− (r, η)) K e D′ . K 1 + τe− (r, η) Moreover, from (16) we also have eD − K e D′ = σ 2 (KD − KD′ ) = σ 2 K



c(x′ ).c(x′ )T c(x).c(x)T − v(x′ ) + r2 v(x) + r2

 ,

eD − K e D′ ) ≤ 2. Thus, exactly as in the determinant argument from the proof of Proposition 1, and hence rank(K we get eα) det(K 1 log − ≤ 2 ψα (e τ− (r, η)) , e D )1−α det(K e D ′ )α 2(α − 1) det(K e α := αK e D′ + (1 − α)K e D . As before, in the singular case, this determinant expression is understood on where K the common support, with determinants replaced by pseudodeterminants. Since V− (r) ≤ Vn (r), and the map V 7→

σ2 V σ 2 r2 + η 2 (v(x) + r2 )

is non-decreasing on V ∈ [0, ∞), we have τe− (r, η) ≤ τen (r, η). Using that ψα is non-decreasing in its argument, the determinant term is bounded by 2 ψα (e τn (r, η)). It remains to bound the mean term. By Lemma A.1, Kα = αKD′ + (1 − α)KD ⪰

r2 − (α − 1)V− (r) C. V− (r) + r2

e α = σ 2 Kα + η 2 KT = σ 2 Kα + η 2 (C + H) satisfies Therefore K e α ⪰ γ C, K

γ := σ 2

r2 − (α − 1)V− (r) + η2 V− (r) + r2

(17)

Tomasz Maciazek

The assumed range of 1 < α < 1 + 1/e τn (r, η) ensures that the scalar prefactor γ > 0. Since with η ̸= 0 the matrix e α may no longer have the same nullspace as C, we cannot apply the Moore-Penrose pseudoinverse directly to K (17) as we did in the proof of Proposition 1. However, for any matrix A and any u ∈ range(A) we have  uT A† u = sup 2uT v − v T Av . v

e α ) (by the bound (17)), we can take A = K e α , C and u = δT to Since δT ∈ range(C) also belongs to the range(K get n o  e α† δT = sup 2δTT v − v T K e α v ≤ sup 2δTT v − v T (γC) v = γ −1 δTT C † δT . δTT K v

v

Taking the supremum over neighbouring pairs and finite test sets, using V− (r) ≤ Vn (r), and plugging the expression for γ gives the following final bound for the Gaussian Rényi mean term  α Vn (r) + r2 α T e† δ K δT ≤ ∆2 (r). 2 T α 2 [σ 2 (r2 − (α − 1)Vn (r)) + η 2 (Vn (r) + r2 )] n   e D (XT )∥Π e D′ (XT ) proves the claim. Combining the determinant and mean bounds for Dα Π

B

BOUNDS FOR Vn AND ∆n

We first present the general bounds for Vn . Proposition 3 (Upper bound for Vn (r).). Let r > 0 and let k : ΩX × ΩX → R be a positive semidefinite kernel satisfying supx∈ΩX k(x, x) = 1. Suppose that k admits a nonnegative lower bound κ := inf x,x′ ∈ΩX k(x, x′ ) ≥ 0. Then, for every n ≥ 1 n−1 Vn (r) ≤ V̄n (r) := 1 − κ2 . n − 1 + r2 If, additionally k(x, x) = 1 for all x ∈ ΩX , then this upper bound is exact, namely Vn (r) = V̄n (r). Proof. By the assumptions, we have 0 ≤ k(x, x) ≤ 1 for all x ∈ ΩX . Define an :=

n−1 . n − 1 + r2

The case n = 1 is immediate. Since the common core D− is empty, V1 (r) = supx∈ΩX k(x, x) = 1 = 1 − a1 κ2 , because a1 = 0. Next, consider n > 1 and denote n− := n − 1. Then D− is of size n− and denote it’s covariate set as X− = {x1 , . . . , xn− }. For x ∈ ΩX , write K := K(X− , X− ),

kx := K(X− , x),

A := K + r2 I.

Then, kD− (x, x′ ) = k(x, x′ ) − kTx A−1 kx′ . We first prove the upper bound. Let 1 ∈ Rn− denote the all-ones vector. Since A is positive definite, we can apply Cauchy-Schwarz inequality with respect to the A-inner product ⟨u, v⟩A = uT Av with u = 1 and v = A−1 kx which gives (1T kx )2 ≤ (1T A1)(kTx A−1 kx ), and thus kTx A−1 kx ≥

(1T kx )2 (1T kx )2 = T . T 1 A1 1 K1 + r2 n−

By the kernel lower bound, 1T kx =

n− X

k(xi , x) ≥ n− κ.

i=1

Since κ ≥ 0, this implies (1T kx )2 ≥ n2− κ2 , thus kTx A−1 kx ≥ κ2

n2− . 1T K1 + r2 n−

Differential Privacy of Gaussian Process Posterior Sampling

Next, since k is positive semidefinite and normalised, we have |k(xi , xj )| ≤ 1T K1 =

n− X

p k(xi , xi )k(xj , xj ) ≤ 1, and hence

k(xi , xj ) ≤ n2− .

i,j=1

It follows that kTx A−1 kx ≥

n2− κ2 n− κ2 = = an κ2 . n2− + r2 n− n− + r2

Substituting into the expression for kD− (x, x) gives kD− (x, x) ≤ k(x, x) − an κ2 ≤ 1 − an κ2 . Taking the supremum over x and over all common cores gives Vn (r) ≤ 1 − an κ2 . Now assume that k(x, x) =  1 for all x ∈ ΩX . We prove the reverse inequality. By definition of κ, there exists a (j) (j) sequence of pairs x1 , x2 such that   j→∞ (j) (j) −−−→ κ. k x1 , x2 (j)

(j)

For each j, let the common core D− consist of n− copies of x1 . Then   (j) (j) K X− , X− = 11T ,

  (j) kx = k x, x1 1,

A = 11T + r2 I.

By the Sherman–Morrison formula, 11T + r2 I

−1

=

1 1 I− 2 11T . 2 r r (n− + r2 )

Consequently, for any x, x′ ∈ ΩX ,         −1 (j) (j) (j) (j) (j) kD− (x, x′ ) = k(x, x′ ) − k x, x1 k x′ , x1 1T 11T + r2 I 1 = k(x, x′ ) − an k x, x1 k x′ , x1 .   (j) (j) (j) Taking x = x′ = x2 , and using k x2 , x2 = 1, gives    2 j→∞ (j) (j) (j) (j) (j) kD− x2 , x2 = 1 − an k x1 , x2 −−−→ 1 − an κ2 . Thus, Vn (r) ≥ 1 − an κ2 . Combining this with the upper bound proves Vn (r) = 1 − an κ2 . The problem of finding bounds for ∆n (r) is more subtle. Below, we treat several special cases which illustrate how ∆n (r) can exhibit a range of different behaviours with respect to n. Throughout, we restrict our considerations to response models with the response bounded in absolute value by a constant MY > 0. We start with the generic upper bound that holds for general bounded-response models. Lemma B.1 (Generic bounded-response bound for ∆n (r)). Assume k(x, x) ≤ 1 for all x ∈ ΩX , and suppose that all training responses satisfy |yi | ≤ MY . Define  1   2, Vn (r) ≥ r2 ,  u 4r Φn (r) := sup = 2 2 Vn (r)  2 0≤u≤Vn (r) (u + r )   2 , Vn (r) ≤ r . (Vn (r) + r2 ) Then

 ∆n (r) ≤ 2MY

1+

n−1 r



p

Φn (r).

Tomasz Maciazek

Proof. Fix neighbouring datasets D = D− ∪ {(x, y)}, D′ = D− ∪ {(x′ , y ′ )}. Let X− and y− denote the commoncore covariates and responses, and write A := K(X− , X− ) + r2 I. For any x0 ∈ ΩX , the common-core posterior mean reads µD− (x0 ) = K(x0 , X− )A−1 y− . By Cauchy–Schwarz in the A−1 -inner product, |µD− (x0 )| ≤

p

q T A−1 y . K(x0 , X− )A−1 K(X− , x0 ) y− −

(18)

Since kD− (x0 , x0 ) ≥ 0, we have K(x0 , X− )A−1 K(X− , x0 ) ≤ k(x0 , x0 ) ≤ 1. Additionally, since A ⪰ r2 I and by the response boundedness we also have T −1 y− A y− ≤

Plugging this into (18) we obtain |µD− (z)| ≤ MY

(n − 1)MY2 ∥y− ∥22 ≤ . 2 r r2 n − 1/r. Thus,

|y − µD− (x)|, |y ′ − µD− (x′ )| ≤ MY

 1+

n−1 r

 .

By the one-step update formula (9) we have µD − µD′ = η kD− (·, x) − η ′ kD− (·, x′ ), for suitable scalars η, η ′ . Putting kD− (x, x) = v and kD− (x′ , x′ ) = v ′ the RKHS norm can be upper bounded as √ √ ∥µD − µD′ ∥H− ≤ ∥ηkD− (·, x)∥H− + ∥η ′ kD− (·, x′ )∥H− = |η| v + |η ′ | v ′ √ √ v v′ ′ ′ + |y − µD− (x )| , = |y − µD− (x)| 2 ′ v(x) + r v(x ) + r2 where in the last step we have applied the formulae for η, η ′ from (9). Since 0 ≤ v, v ′ ≤ Vn (r), we have √ √ p v v′ , ≤ Φn (r). 2 ′ 2 v(x) + r v(x ) + r Thus ∥µD − µD′ ∥H− ≤ 2MY

√   n−1 p Φn (r). 1+ r

Taking the supremum over neighbouring pairs gives the desired result. The Lemma below shows that the

n-growth of the worst-case upper bound is in fact tight.

Lemma B.2. Let ΩX = [0, 1] and suppose that all training responses satisfy |yi | ≤ MY . There exists a positive semidefinite kernel k : ΩX × ΩX → R satisfying supx∈ΩX k(x, x) ≤ 1 for which ∆n (r) has the following lower bound.  √ MY 1 + r 2 + n − 1 √ ∆n (r) ≥ . r 1 + r2 (2 + r2 ) Proof. Let (ei )i≥1 denote the standard basis of ℓ2 . Choose mutually distinct points x0 ,

a1 , a2 , . . . ,

b1 , b2 , . . .

in [0, 1]. Define a feature map φ : [0, 1] → ℓ2 by m

φ(x0 ) := 0,

φ(ai ) := ei ,

1 X φ(bm ) := √ ei , m i=1

and set φ(x) := 0 for all remaining x ∈ [0, 1]. Now define the GP kernel through this feature map k(x, x′ ) := ⟨φ(x), φ(x′ )⟩ℓ2 . Then k has bounded diagonal since k(x, x) = ∥φ(x)∥2ℓ2 ≤ 1 for every x ∈ [0, 1].

Differential Privacy of Gaussian Process Posterior Sampling

Fix n ≥ 2. Consider the common-core dataset D− := {(ai , MY )}n−1 i=1 . Then K(X− , X− ) = I and K(X− , bn−1 ) = √ 1 1. Thus, the common-core posterior mean at x = bn−1 is n−1 2

µD− (bn−1 ) = K(bn−1 , X− ) K(X− , X− ) + r I

−1

√ MY n − 1 y− = . 1 + r2

The corresponding common-core posterior variance is v := kD− (bn−1 , bn−1 ) = k(bn−1 , bn−1 ) − K(bn−1 , X− ) K(X− , X− ) + r2 I

−1

K(X− , bn−1 )

2

=1−

r 1 = . 1 + r2 1 + r2

Now define neighbouring datasets D := D− ∪ {(bd , −MY )}, D′ := D− ∪ {(x0 , MY )}. Since φ(x0 ) = 0, we have k(·, x0 ) ≡ 0 and hence kD− (·, x0 ) ≡ 0. Thus, by the one-step update formula (9) µD − µD′ = η kD− (·, bn−1 ),

η=

−MY − µD− (bn−1 ) . v(x) + r2

Using the expressions above, we have √ √ −1   2 n−1 MY (1 + r2 + n − 1) r 2 η = −MY 1 + +r =− . 1 + r2 1 + r2 r2 (2 + r2 ) By the RKHS property we have we have ∥µD − µD′ ∥2H− = η 2 kD− (bn−1 , bn−1 ) = η 2 v. Substituting the expressions for η and v gives  √ MY 1 + r 2 + n − 1 √ . ∥µD − µD′ ∥H− = r 1 + r2 (2 + r2 ) The claimed lower bound follows by taking the supremum over neighbouring datasets D ∼ D′ . Lemma B.3. Let k : ΩX × ΩX → R be a positive semidefinite kernel satisfying supx∈ΩX k(x, x) ≤ 1. Assume that there exist admissible datapoints (x0 , y0 ) and (x1 , y1 ) such that y0 > 0, y1 < 0 and k(x0 , x0 ) = 1,

k(x0 , x1 ) > 0,

k(x1 , x1 ) − k(x0 , x1 )2 > 0.

Then p |y1 | k(x1 , x1 ) − k(x0 , x1 )2 ∆n (r) ≥ > 0. 1 + r2 Proof. Let D− consist of n − 1 repeated copies of (x0 , y0 ) and set D = D− ∪ {(x0 , y0 )}, D′ = D− ∪ {(x1 , y1 )} and define n n−1 , γn := an := 2 n−1+r n + r2 and let XT = {x1 }. Let D− consist of n − 1 repeated copies of (x0 , y0 ). Then, K(X− , X− ) = 11T and using the Sherman-Morrison formula it is straightforward to verify that µD− (x1 ) = an k(x0 , x1 ) y0 ,

µD (x1 ) = γn k(x0 , x1 ) y0 ,

vn := kD− (x1 , x1 ) = k(x1 , x1 ) − an k(x0 , x1 )2 > 0.

Also, because an ≥ 0 and supx k(x, x) ≤ 1, we have vn ≤ k(x1 , x1 ) ≤ 1. For D′ , the one-point update formula (9) gives µD′ (x1 ) = µD− (x1 ) + kD− (x1 , x1 )

y1 − µD− (x1 ) y1 − an k(x0 , x1 ) y0 = an k(x0 , x1 ) y0 + vn kD− (x1 , x1 ) + r2 vn + r 2 =

an k(x0 , x1 ) y0 r2 + vn y1 . vn + r 2

Tomasz Maciazek

Hence, an k(x0 , x1 ) y0 r2 + vn y1 vn + r 2 vn (γn k(x0 , x1 ) y0 − y1 ) + k(x0 , x1 ) y0 r2 (γn − an ) . = vn + r 2

µD (x1 ) − µD′ (x1 ) = γn k(x0 , x1 ) y0 −

Since k(x0 , x1 ) > 0, y0 > 0, y1 < 0, and γn − a n =

r2 (n + r2 )(n − 1 + r2 )

≥ 0,

both terms in the numerator are nonnegative. Thus µD (x1 ) − µD′ (x1 ) ≥

vn |y1 | vn (γn k(x0 , x1 ) y0 − y1 ) ≥ . vn + r 2 vn + r 2

By the reproducing property in H− and Cauchy-Schwarz |µD (x1 ) − µD′ (x1 )| = ⟨µD − µD′ , kD− (x1 , ·)⟩H− ≤ ∥µD − µD′ ∥H− Therefore ∥µD − µD′ ∥2H− ≥

µD (x1 ) − µD′ (x1 ) vn

q

√ kD− (x1 , x1 ) = ∥µD − µD′ ∥H− vn .

2 ≥ |y1 |2

vn . (vn + r2 )2

Since k(x1 , x1 ) − k(x0 , x1 )2 ≤ vn ≤ 1, we have vn k(x1 , x1 ) − k(x0 , x1 )2 ≥ . (vn + r2 )2 (1 + r2 )2 Hence

p |y1 | k(x1 , x1 ) − k(x0 , x1 )2 ∥µD − µD′ ∥H− ≥ . 1 + r2 Taking the supremum over admissible neighbouring pairs gives proves the claim. On the other hand, there exist typical situations where ∆n (r) is upper-bounded by a constant or even decreases with n. Lemma B.4 (Constant kernel). Let k(x, x′ ) ≡ 1 for all x, x′ ∈ ΩX and suppose that all admissible responses satisfy |y| ≤ MY . Then √ √ 2MY r2 + n − 1 ∆n (r) ≤ = O(1/ n). r(r2 + n) Proof. Let D = D− ∪ {(x, y)}, D′ = D− ∪ {(x′ , y ′ )} where D− is any common core of size n − 1. Since the kernel is constant for any X we have K(X, X) = 11T and it is straightforward to verify (using Sherman–Morrison formula) that for any dataset D of size n the posterior mean is a constant function and that the posterior covariance is also constant and given by Pn yi r2 kD− (u, v) = 2 . µD (·) = 2i=1 , r +n r +n−1 Thus, µD − µD′ = (y − y ′ )/(r2 + n) and H− consists of constant functions. A constant function g has RKHS norm √ r2 + n − 1 ∥g∥H− = |g| . r Hence √ |y − y ′ | r2 + n − 1 ∥µD − µD′ ∥H− = 2 . r +n r Taking the supremum over |y|, |y ′ | ≤ MY gives |y − y ′ | ≤ 2MY from which the final result follows.

Differential Privacy of Gaussian Process Posterior Sampling

Lemma B.5 (Purely diagonal kernel). Let k(x, x′ ) = q(x) if x = x′ with 0 ≤ q(x) ≤ 1 and zero otherwise. Suppose that all admissible responses satisfy |yi | ≤ MY . Then √ 2MY    , 0 < r ≤ 1, √r ∆n (r) ≤    2 2MY , r ≥ 1. 1 + r2 Proof. Let D− be the common core. For u ∈ ΩX , define the overlap index set I∩ (u) := {i : xi = u, (xi , yi ) ∈ D− },

n∩ (u) := #I∩ (u),

S∩ (u) :=

X

yi .

i∈I∩ (u)

For the purely diagonal kernel, the common-core posterior kernel is again diagonal. Its diagonal value is v(u) := kD− (u, u) =

q(u)r2 ≤ q(u) ≤ 1. ∩ (u)q(u)

r2 + n

The common-core posterior mean is bounded as |µD− (u)| =

|S∩ (u)| q(u)|S∩ (u)| ≤ ≤ MY . n∩ (u) ∩ (u)q(u)

r2 + n

where we have used the fact that |S∩ (u)| ≤ n∩ (u)MY . Hence, for any additional admissible observation (u, y) we have |y − µD− (u)| ≤ 2MY . Now let D = D− ∪ {(x, y)} and D′ = D− ∪ {(x′ , y ′ )}. If x ̸= x′ , then kD− (x, ·) and kD− (x′ , ·) are orthogonal in H− , because kD− is diagonal. Then, the one-point update formula (9) yields √ ∥µD − µD′ ∥H− ≤ 2 2MY sup

0≤v≤1 r

v

2+v

.

It is straightforward to check that the same bound holds also when x = x′ . The final result is obtained from the fact that  1  √  , r2 ≤ 1, v 2r sup 2 = 1  0≤v≤1 r + v  2 , r2 > 1. r +1 Lemma B.6. Let ΩX ⊆ R and let k(x, x′ ) = exp (−|x − x′ |/ℓ) with ℓ > 0. Suppose that all admissible responses satisfy |yi | ≤ MY . Then, for every admissible dataset D and for all x ∈ R we have |µD (x)| ≤ MY . Proof. For D = (X, y) write X = (x1 , . . . , xn ) and assume that x1 < x2 < · · · < xn . Denoting K := K(X, X) we have µD (X) = K(K + r2 I)−1 y = KK −1 (I + r2 K −1 )−1 y = (I + r2 K −1 )−1 y. For i < j we can represent the entries of K as  Kij = exp

xi − x j ℓ

 =

j−1 Y k=1

 ρk ,

ρk := exp

xk − xk+1 ℓ

 .

One can verify by a straightforward calculation that K −1 is tridiagonal with entries given by  −1  K = 11

1 , 1 − ρ21

 −1  1 1 ρ2i , K = + , 2 2 ii 1 − ρn−1 1 − ρi−1 1 − ρ2i  −1    ρi K = K −1 i+1,i = − . i,i+1 1 − ρ2i

 −1  K = nn

2 ≤ i ≤ n − 1,

Tomasz Maciazek

Using these formulas we also get that K −1 has nonpositive off-diagonal entries and non-negative row-sums. Consequently, A := I + r2 K −1 has nonpositive off-diagonal entries and row-sums ≥ 1. Now find index i∗ such that µD (xi∗ ) = maxi µD (xi ). Then, for any j we have µD (xj ) ≤ µD (xi∗ ) Since Aij ≤ 0 for i ̸= j, we also have Ai∗ ,j µD (xj ) ≥ Ai∗ ,j µD (xi∗ ). Summing over j yields yi∗ = [A µD (X)]i∗ =

X

Ai∗ ,j µD (xj ) ≥ µD (xi∗ )

j

X

Ai∗ ,j ≥ µD (xi∗ ),

j

P since the row-sum satisfies j Ai∗ ,j ≥ 1. By the response-boundedness, this yields µD (xj ) ≤ µD (xi∗ ) ≤ MY . By picking i∗ such that −µD (xi∗ ) = maxi (−µD (xi )) and repeating the reasoning above, we get −µD (xj ) ≤ −µD (xi∗ ) ≤ MY , i.e. µD (xj ) ≥ −MY and consequently |µD (xj )| ≤ MY

for all

xj ∈ X.

It remains to show that |µD (x)| ≤ MY for any x outside the training set X. Assume first that xi ≤ x ≤ xi+1 for some 1 ≤ i < n. We claim that |µD (x)| is convex on the interval [xi , xi+1 ], thus max

xi ≤x≤xi+1

|µD (x)| = max{|µD (xi )|, |µD (xi+1 )|} ≤ MY .

To see this, note that µD (x) as a function of x is a linear combination of e−x/ℓ and ex/ℓ , so it satisfies the equation ∂x2 µD (x) = µD (x)/ℓ2 . In the region where µD (x) ≥ 0 this gives ∂x2 |µD (x)| > 0 i.e. |µD (x)| is convex. In the region where µD (x) < 0 we have |µD (x)| = −µD (x) and thus ∂x2 |µD (x)| = −µD (x)/ℓ2 > 0 as well. The same conclusion holds on the exterior intervals (if they are present). Indeed, on the left exterior interval ΩX ∩ (−∞, x1 ), the function µD (x) is of the form µD (x) = a− ex/ℓ for some constant a− , and hence |µD (x)| is maximised at x1 if x1 is finite, while it tends to zero as x1 → −∞ when ΩX is unbounded to the left. Analogous reasoning holds for the right exterior interval. When the dataset D contains repeated covariates, then we consider the collapsed dataset D̄ = (X̄, ȳ) of the unique covariates x̄1 < x̄2 < · · · < x̄m where x̄k appears nk times in the original D, nk ≥ 1. Suppose Pnk that yk,1 , . . . , yk,nk are the responses of x̄k in D. Then, we define the collapsed response of x̄k as ȳk := n1k i=1 yk,i . In this notation, at the unique training points we have µD (X̄) = µ̄D̄ (X̄) := (I + RK̄ −1 )−1 ȳ, where R is a diagonal matrix with entries Rkk = r2 /nk and K̄ = K(X̄, X̄). The rest of the proof follows exactly the same as before with µD replaced with µ̄D . Lemma B.7 (Exponential kernel in 1D). Let ΩX ⊆ R be a possibly unbounded interval and let k(x, x′ ) = exp (−|x − x′ |/ℓ), x, x′ ∈ ΩX with ℓ > 0. Suppose that all admissible responses satisfy |yi | ≤ MY . Then, for every r > 0,  1   Vn (r) ≥ r2 ,  4r2 , p u ∆n (r) ≤ 4MY Φn (r), Φn (r) = sup = 2 2 Vn (r)  2 0≤u≤Vn (r) (u + r )   2 , Vn (r) ≤ r . (Vn (r) + r2 ) Proof. Let D− = (X− , y − ) be the common core of D ∼ D′ . Denote by µD− and kD− the corresponding posterior mean and covariance and set v(x) = kD− (x, x). By the one-step update formula 9 we have µD− ∪{(x,y)} − µD− =

y − µD− (x) kD− (x, ·). r2 + v(x)

p The reproducing property gives v(x) = ∥kD− (x, ·)∥H− . The 1D-exponential kernel has the property that |µD (x)| ≤ MY for any D and any x, see Lemma B.6. Hence, y − µD− (x) 2MY ≤ 2 . 2 r + v(x) r + v(x)

Differential Privacy of Gaussian Process Posterior Sampling

Set D = D− ∪ {(x, y)} and D′ = D− ∪ {(x′ , y ′ )}. By the one-step update formula 9 we have   √ p ∥kD− (x′ , ·)∥H− ∥kD− (x, ·)∥H− v + ≤ 4M sup = 4MY Φn (r). ∥µD − µD′ ∥H− ≤ 2MY Y 2 2 ′ 2 r + v(x) r + v(x ) 0≤v≤Vn (r) r + v Taking the supremum over D ∼ D′ yields the result. Lemma B.8. Let k be a positive semidefinite normalised kernel with bounded diagonal with RKHS Hk and let r > 0. Suppose that the responses are generated by a fixed function f∗ ∈ Hk i.e., for every admissible dataset D = (X, y) with X = (x1 , . . . , xn ), we have yi = f∗ (xi ), i = 1, . . . , n. Let µD denote the GP posterior mean as defined in Section 2. Then Vn (r) ∆n (r) ≤ 2∥f∗ ∥Hk 2 . r + Vn (r) Proof. Let D = D− ∪ {(x0 , f∗ (x0 ))}, D′ = D− ∪ {(x′0 , f∗ (x′0 ))}. Write X− for the covariates in D− . Define the common-core posterior kernel as in Section 2  −1 (r) kD− (x, x′ ) := k(x, x′ ) − kX− (x)T KX− X− kX− (x′ ),

(19)

(r)

with KX− X− = K(X− , X− ) + r2 I, kX− (x) = K(X− , x). Also denote v(x) := kD− (x, x). First note that kD− (x, ·) ∈ Hk since from the representation (19) we have that kD− (x, ·) is a finite linear combination of k(x′ , ·) with x′ ∈ X− . By the one-point update formula (9), we have µD − µD′ =

f∗ (x′0 ) − µD− (x′0 ) f∗ (x0 ) − µD− (x0 ) k (x , ·) − kD− (x′0 , ·) ∈ H− . D 0 − r2 + v(x0 ) r2 + v(x′0 )

(20)

We now bound each term. Since kD− (x, ·) ∈ Hk , we have  f∗ , kD− (x, ·) H = k

 −1 (r) f∗ , k(x, ·) − kX− (x)T KX− X− kX− (·)

= f∗ (x) − kX− (x)T



(r) KX− X−

−1

Hk

f∗ (X− ) = f∗ (x) − µD− (x),

where we have used the RKHS property ⟨f∗ , k(x, ·)⟩Hk = f (x). Thus, by Cauchy–Schwarz, |f∗ (u) − µD− (u)| ≤ ∥f∗ ∥Hk ∥kD− (u, ·)∥Hk .  −1 (r) Next, denoting a = KX− X− kX− (x), we get (r)

∥kD− (x, ·)∥2Hk = k(x, x) − 2aT kX− (x) + aT K(X− , X− )a = v(x) − aT kX− (x) + aT KX− X− a − r2 ∥a∥2 = v(x) − r2 ∥a∥22 ≤ v(x). Thus, for any x we have |f∗ (x) − µD− (x)| ≤ ∥f∗ ∥Hk

p v(x) which in turn yields

p f∗ (x) − µD− (x) ∥f∗ ∥Hk v(x) p v(x) v(x) = ∥f∗ ∥Hk 2 , ∥kD− (x, ·)∥H− ≤ r2 + v(x) r2 + v(x) r + v(x) p where we have also used the reproducing property ∥kD− (x, ·)∥H− = v(x). Applying this estimate twice to (20) with x = x0 and x = x′0 respectively and using the triangle inequality in H− gives     v(x0 ) v(x′0 ) v Vn (r) ∥µD − µD′ ∥H− ≤ ∥f∗ ∥Hk + ≤ 2∥f ∥ sup = 2∥f∗ ∥Hk 2 . ∗ Hk 2+v r2 + v(x0 ) r2 + v(x′0 ) r r + Vn (r) 0≤v≤Vn (r) Taking the supremum over neighbouring datasets D ∼ D′ gives the final result.

Tomasz Maciazek

C

FURTHER DETAILS OF THE MEMBERSHIP ATTACK

We evaluate membership leakage using a LiRA-style shadow-model attack. The attack targets a fixed point x0 and distinguishes between two hypotheses: Hin : x0 ∈ D,

Hout : x0 ∈ / D.

For each hypothesis, the adversary generates shadow datasets from the same data-generating distribution used in the experiment (uniform on ΩX = [0, 1]). In the step-function experiment, responses are generated from ( −1, x < 1/2, fstep (x) = 1, x ≥ 1/2, with bounded additive label noise yi = (1 − Mξ )fstep (xi ) + ξi ,

ξi ∼ Unif[−Mξ , Mξ ].

Under Hin , the target point x0 is inserted into the dataset and the remaining points are sampled from the uniform distribution on [0, 1]. Under Hout , all n points are sampled from the uniform distribution on [0, 1]. For a shadow dataset D = (X, y) the posterior mean and normalized posterior variance at x0 are µD (x0 ) = k(x0 , X) K(X, X) + r2 I and vD (x0 ) = 1 − k(x0 , X) K(X, X) + r2 I ′

−1

−1

y,

k(X, x0 ),

where k(x, x ) = exp(−|x − x |/ℓ) is the exponential kernel function. The released object consists of L posterior sample-path values at the target point, p (ℓ) fD (x0 ) = µD (x0 ) + σ vD (x0 )Zl , Zl ∼ N (0, 1), l = 1, . . . , L. (21) The attack uses two scalar statistics based on L released values: 1 fbD (x0 ) = L

L X

L

(ℓ)

fD (x0 ),

vbD (x0 ) =

ℓ=1

2 1 X  (ℓ) fD (x0 ) − fbD (x0 ) . 2 Lσ ℓ=1

Thus fbD (x0 ) captures the posterior mean signal available from the released paths, while vbD (x0 ) captures the posterior covariance signal. For L = 1, the empirical variance statistic is undefined and the attack uses only fbD (x0 ). The LiRA score is obtained by fitting separate density models for the in and out shadow distributions. Let (f ) (f ) ρin , ρout denote the fitted densities of fbD (x0 ) under Hin and Hout respectively. The sample-mean statistic is modelled by a two-component latent Gaussian mixture. Specifically, the model introduces a latent Gaussian mixture variable ϕ∼

2 X

ωj N (mj , τj2 ),

 N (mj , τj2 ) : ϕ 7→ ρGauss ϕ; mj , τj2 :=

j=1

1 √

τj 2π

  exp −(ϕ − mj )2 /(2τj2 )

and models the observed statistic as fbD (x0 ) = u + η,

u = tanh(ϕ/2),

η ∼ N (0, s2 ).

This latent model is motivated by Equation 21 wherepµD (x0 ) takes values in [−1, 1] (by Lemma B.6) and thus is modelled by tanh(ϕ/2) whereas the contribution σ vD (x0 )Zl is modelled by the random normal variable η. Consequently, the fitted density is obtained by marginalising over the distribution of the latent variable ϕ ρ

(f )

Z ∞ ρ(t | ϕ)ρ(ϕ)dϕ =

(t) = −∞

2 X j=1

Z ∞ ωj −∞

  ρGauss t; tanh(ϕ/2), s2 ρGauss ϕ; mj , τj2 dϕ,

Differential Privacy of Gaussian Process Posterior Sampling

Figure 5: MIA success and (ε, δ)-DP bounds (n = 10, δ = 0.05). Left: Excess TPR@FPR = 0.01 vs. ε for different σ (L = 1). Error bands show standard deviation over 10 random seeds. Right: Effect of releasing L posterior sample paths. The LiRA attack starts with random-guess effectiveness (T P R ≈ F P R) at L = 1 and strengthens with L approaching the non-private benchmark while the (ε, δ)-DP bound grows sub-linearly with L. Parameters: r = 1, σ = 5, n = 10, δ = 0.05. Error bands show standard deviation over 10 random seeds. with the parameters fitted separately for the in- and out- shadow distributions by maximum likelihood. The integral is evaluated numerically. For the nonnegative sample-variance statistic vbD (x0 ) the attack instead fits a two-component Gaussian mixture after the logarithmic transformation i.e., ϕv = log vbD (x0 ),

ϕv ∼

2 X

 ω̃j N m̃j , τ̃j2 .

j=1

Given a fresh evaluation release, the fitted likelihood-ratio scores are     (f ) (f ) Sf = log ρin fbD (x0 ) − log ρout fbD (x0 ) , and, for L > 1, (v)

(v)

Sv = log ρin (b vD (x0 )) − log ρout (b vD (x0 )) . The combined attack score is S = Sf + S v . Equivalently, this is a plug-in product approximation to the joint likelihood ratio of (fbD (x0 ), vbD (x0 )). When L = 1 the combined score reduces to S = Sf . Figure 6 shows that the latent models model the empirical in- and out- histograms well. Attack performance is evaluated on independent in- and out- evaluation datasets. We report the ROC curve and true-positive rates at fixed false-positive rates. For a score threshold t the Neyman–Pearson test declares membership when S ≥ t. This gives false-positive and true-positive rates FPR(t) = Pout {S ≥ t},

TPR(t) = Pin {S ≥ t}.

To report TPR@α we vary t until the false-positive rate under the out distribution reaches the target level α, and then evaluate the true-positive rate under the in distribution at the same threshold. Thus, for a threshold tα satisfying FPR(tα ) = α we report TPR@α = TPR(tα ) = Pin {S ≥ tα }. This evaluation measures how well an adversary equipped with shadow samples from the two membership hypotheses and knowledge of the release parameters (ℓ, r, σ, L) can distinguish whether the target point was included in the training dataset. Figure 5 presents results for TPR@FPR = 0.01.

D

FURTHER DETAILS OF THE PRIVACY-UTILITY EXPERIMENTS

Details of the 1D excursion-set utility experiment. For each experimental configuration we generate a   (j) 3 (j) pool of Npair = 10 independent pairs f∗ , D , j = 1, . . . , Npair on ΩX = [0, 1]. For pair j the covariates

Tomasz Maciazek

Figure 6: Shadow-distribution histograms for the LiRA-style attack at x0 = 1/2 showing the sample-path mean statistic fbD (x0 ) in the top row and the log-variance statistic log vbD (x0 ) in the bottom row. Each column corresponds to a different choice of the effective ridge parameter r and kernel scale σ with the corresponding approximate DP value ε shown in the figure title. The comparison illustrates that privacy depends sharply on effective ridge regularisation. When r is small, the in- and out-shadow distributions remain well separated, and even a very large posterior sampling scale σ does not yield a small privacy parameter. In contrast, increasing r substantially reduces the distinguishability of the two hypotheses: in the right column the in/out histograms and fitted densities nearly overlap, corresponding to much stronger membership privacy protection. Dataset size |D| = 10, δ = 0.05, the LiRA attack uses L = 3 independent posterior draws. The responses are generated as noiseless (Mξ = 0). (j)

(j)

x1 , . . . , xn are sampled i.i.d. from Unif[0, 1]. The latent target is sampled from a centred GP with exponential kernel e k(x, x′ ) = exp(−|x − x′ |/ℓgen ) with ℓgen = 1. The sampled path is rescaled to have supremum norm 1 − Mξ and the observations for pair j are generated as   (j) (j) (j) (j) (j) i.i.d. xi + ξi , ξi ∼ Unif[−Mξ , Mξ ]. yi = f∗ (j)

(j)

More precisely, if the unnormalised prior GP draw is denoted by f˜∗ , then f∗

(j) (j) = (1 − Mξ )f˜∗ / fe∗

. Thus

the amplitude noise-to-signal ratio is N SR = Mξ /(1 − Mξ ). In all reported experiments we take the excursion threshold t = 0, so the target excursion set is o   n (j) (j) Ωt f ∗ = x ∈ ΩX : f∗ (x) ≥ t . grid Norm ∥ · ∥∞ is approximated on a fixed uniform evaluation grid xgrid with m = 800 and all integrals 1 , . . . , xm over ΩX are approximated using the trapezoidal rule on the same grid. We write

s∗ (x) := 1 {f∗ (x) ≥ t} . where 1{·} is the indicator function.   (j) For a candidate hyperparameter triple θ = (ℓ, r, σ), we fit the GP posterior for each pair f∗ , D(j) and compute

Differential Privacy of Gaussian Process Posterior Sampling

the posterior excursion probability µD (x) − t p pD,θ (x) := Pf ∼ΠD,θ [f (x) ≥ t] = Φ σ kD (x, x)

! ,

where Φ is the CDF of the univariate normal distribution. The hyperparameter-search objective is the prioraverage integrated BCE LbBCE (θ) =

1

Npair

X

Npair j=1

Z

  (j) LBCE ΠD(j) ,θ , f∗ ,

LBCE (ΠD,θ , f∗ ) :=

ℓBCE (s∗ (x), pD (x)) dx. ΩX

The unconstrained benchmark hyperparameters are selected by θ∗ ∈ arg min LbBCE (θ), θ

whereas the privacy-feasible hyperparameters are selected by the constrained search priv θL ∈ arg

min θ: εL (θ)<ε0

LbBCE (θ).

Thus the constrained search differs from the unconstrained search only through the privacy feasibility constraint. Both searches optimise the same BCE objective. The DP bound εL (θ) is evaluated using the exponential-kernel sensitivity formula from Lemma B.7 together with Theorem 3.1 and the RDP composition rule (3) numerically optimised over α with L being the number of released posterior sample paths. The hyperparamater search is performed on a coarse grid over (ℓ, r, σ) followed by two rounds of grid refinement around the top-3 best average-BCE candidates with the private refinement restricted to feasible candidates satisfying εL (θ) < ε0 . priv , we evaluate both the posterior probability map and the actual randomized release. After selecting θ∗ and θL The deterministic estimate associated with a candidate θ is

sbθ (x) := 1 {pD,θ (x) ≥ C} , where the cutoff C is selected on a separate validation set by maximising expected IoU where R b θ ∩ Ωt (f∗ )  Ω dx s∗ (x)b sθ (x) ΩX b IoU Ωθ , Ωt (f∗ ) = =R . dx 1 {b sθ (x) = 1 or s∗ (x) = 1} b θ ∪ Ωt (f∗ ) Ω ΩX 

  b θ , f∗ is the reported baseline IoU. It is a deterministic post-processing The finite-grid approximation of IoU Ω of the posterior excursion probabilities and is not the randomized posterior-sample release. The utility of the randomized released excursion set is evaluated by Monte Carlo simulation of the actual Lpath mechanism where we draw B independent sets of L posterior sample paths. Namely, for each stored pair   (j) (j) f∗ , D and Monte Carlo repetition b = 1, . . . , B, we draw independent posterior paths (1)

(L) i.i.d.

fj,b , . . . , fj,b

∼ ΠD(j) ,θ

on the evaluation grid. In the experiments we take B = 50. For a fixed collection of posterior sample paths f (1) , . . . , f (L) the indicator function of the released (randomised) excursion set is ( ) L o 1 X n (ℓ) (L) sbθ (x) := 1 1 f (x) ≥ t ≥ c , L ℓ=1

where the vote-fraction cutoff c is also selected on a separate validation set by maximising expected IoU. The b (L) . released excursion set for this collection of posterior sample paths is denoted by Ω θ For fixed pair (f∗ , D), the Monte Carlo mean over the B sample-path set realisations dMC = 1 IoU B (L)

B X b=1

  b (L) , Ωt (f∗ ) IoU Ω b,θ

Tomasz Maciazek

0.1

0.3

NSR 0.11 0.43 ε, unconstrained 2· 104 500 deff , unconstrained 52.7 [52.0, 53.4] 17.7 [17.6, 17.9] deff , ε < 10 2.26 [2.25, 2.26] 2.21 [2.20, 2.22] relative BCE increase 115% [67%, 187%] 33.5% [14.0%, 60.6%] IoU, unconstrained 0.957 [0.932, 0.976] 0.922 [0.877, 0.956] 1-path IoU, ε < 10 0.854 [0.799, 0.895] 0.833 [0.774, 0.880] relative IoU gap 10.1% [7.0%, 14.8%] 8.7% [5.5%, 13.4%]

0.5

0.6

1.00 1.50 63.6 21.4 8.47 [8.42, 8.53] 6.01 [5.97, 6.05] 2.15 [2.14, 2.15] 2.08 [2.07, 2.09] 8.8% [0.2%, 19.4%] 2.2% [0%, 7.4%] 0.874 [0.806, 0.927] 0.844 [0.760, 0.905] 0.797 [0.720, 0.855] 0.772 [0.680, 0.839] 7.7% [4.1%, 13.1%] 7.3% [2.9%, 13.4%]

Table 2: The effect of varying the NSR set size on the released 1D excursion set utility at fixed training set size |D| = 100. Median [interquartile range] over 103 draws of (D, f∗ )-pairs across a range NSR values. The private candidate is the best average-BCE-selected setting satisfying ϵ < 10. (L)

b estimates the expected IoU of the randomized release Ω θ

given (f∗ , D). The reported released-set IoU is   (L) d M C over the stored pairs f∗(j) , D(j) . Similarly, obtained by summarising the median and IQR of the values IoU   (j) the reported L-path IoU SD is the distribution of the MC standard deviations over the stored pairs f∗ , D(j) . Such a standard deviation of the released-set IoU isolates the randomness due to posterior sampling for a fixed   (j) f∗ , D(j) -pair. The BCE increase reported in the final table is  computed  pairwise by comparing the privacy-feasible and uncon(j) (j) strained BCE-selected candidates on the same f∗ , D -pairs:   (j) LBCE ΠD(j) ,θpriv , f∗ L  − 1.  (j) LBCE ΠD(j) ,θ∗ , f∗ The tables report median and interquartile range over the stored pairs for the DP bound, effective dimension deff = tr{K(K +r2 I)−1 }, BCE increase, baseline IoU from the unconstrained hyperparameter search, released-set IoU, and the conditional standard deviation of the released-set IoU. Table 2 shows that as NSR increases the unconstrained utility optimum becomes naturally more regularized: both ε and deff decrease sharply. Consequently, the additional utility cost of the constrained search decreases as NSR grows. The results complement Table 1 from the main text. Dependence of the privacy-utility tradeoff on |D|. Table 3 shows that the cost of the privacy constraint decreases with sample size, as is typical in differential privacy. When |D| is small, each record has high influence and enforcing single-digit privacy requires strong regularisation. In our setting this appears as increased posteriorsample variability. At |D| = 50 imposing ε < 10 reduces deff from 6.63 to 1.13 and yields a relative one-path IoU gap of 8.7%. As |D| grows, the privacy-feasible posterior remains substantially regularised but loses less downstream utility: the relative BCE increase drops from 14.2% at |D| = 50 to 0.6% at |D| = 400, while the relative IoU gap drops from 8.7% to 4.1%. The randomized release also becomes more stable: the one-path IoU standard deviation decreases from 0.082 to 0.038 over the same range. In particular, at |D| = 400 the privacy-feasible posterior reduces deff from 14.6 to 4.25, yet the one-path IoU remains high at 0.874. D.1

Synthetic excursion set experiment in 2D

We consider a synthetic spatial excursion-set problem on ΩX = [0, 1]2 , intended to mimic a smooth pollution field with a small number of broad high-concentration regions. We generate a random pollution field by drawing a small number of anisotropic Gaussian-type plumes with random centres, widths and amplitudes and forming their superposition h(x). More specifically, we draw K anisotropic Gaussian plumes and set K X



 1 ⊤ −1 h(x) = aj exp − (x − xcj ) Σj (x − xcj ) , 2 j=1

2 2 Σj = Rϕj diag(wj,1 , wj,2 )Rϕ⊤j .

Differential Privacy of Gaussian Process Posterior Sampling

|D|

50

100

ε, unconstrained 90.4 deff , unconstrained 6.63 [6.55, 6.70] deff , ε < 10 1.13 [1.12, 1.13] relative BCE increase 14.2% [1.3%, 33.0%] IoU, unconstrained 0.852 [0.790, 0.915] 1-path IoU mean, ε < 10 0.727 [0.629, 0.811] 1-path IoU SD, ε < 10 0.082 [0.055, 0.109] relative IoU gap 8.7% [5.5%, 13.4%]

200

63.6 8.47 [8.42, 8.53] 2.15 [2.14, 2.15] 8.8% [0.2%, 19.4%] 0.874 [0.806, 0.927] 0.797 [0.720, 0.855] 0.056 [0.040, 0.080] 7.7% [4.1%, 13.1%]

400

39.0 16.7 13.0 [12.9, 13.0] 14.6 [14.5, 14.6] 2.17 [2.16, 2.17] 4.25 [4.25, 4.26] 4.3% [0%, 11.5%] 0.6% [0%, 2.7%] 0.888 [0.831, 0.935] 0.913 [0.878, 0.946] 0.835 [0.767, 0.880] 0.874 [0.822, 0.909] 0.051 [0.039, 0.069] 0.038 [0.027, 0.049] 5.7% [2.9%, 10.0%] 4.1% [2.4%, 6.7%]

Table 3: The effect of varying the training set size on the released 1D excursion set utility at fixed N SR = 1. Median [interquartile range] over 103 draws of (D, f∗ )-pairs across a range of training dataset sizes. The private candidate is the best average-BCE-selected setting satisfying ε < 10. The relative IoU gap is computed with respect to the unconstrained IoU. The 1-path IoU mean, SD reported for a sample of 50 posterior draws. Here xcj ∈ [0, 1]2 , aj > 0, wj,1 , wj,2 > 0, ϕj ∈ [0, 2π] are sampled randomly and Rϕ denotes the rotation by angle ϕ. In the presented expetiments we take K = 2. To obtain high- and low-pollution regions, we apply the sigmoid transformation 1 , ge∗ (x) = 1 + exp[−γ(h(x) − bh )] where bh is chosen so that 65% grid values of h are below bh . The resulting field is then normalised to give g∗ : ΩX → [0, 1] on ΩX . Since the GP prior is centred at zero, we finally use the transformed field f∗ (x) = 2g∗ (x) − 1 ∈ [−1, 1] as the latent response function. Given n uniformly sampled sensor locations xi ∈ ΩX , observations are generated as yi = (1 − Mξ )f∗ (xi ) + ξi ,

ξi ∼ Unif[−Mξ , Mξ ].

The target set is the excursion set of the attenuated signal, Ωt = {x ∈ ΩX : (1 − Mξ )f∗ (x) ≥ t},

t = (1 − Mξ )(2q − 1),

where q ∈ (0, 1) is the excursion threshold on the original [0, 1] scale of g∗ . In the presented experiments we take q = 1/2 i.e. t = 0. For each dataset we fit a GP posterior with exponential kernel k(x, x′ ) = exp(−∥x − x′ ∥/ℓ), effective ridge parameter r, and posterior scale σ. As in the 1D experiment, the hyperparameters (ℓ, r, σ) are selected by minimising the mean binary cross-entropy between the true excursion labels and the posterior excursion probabilities ! µD (x) − t p . pD (x) = Φ σ kD (x, x) We perform this selection both without a privacy constraint and under the sufficient privacy constraint ε < ε0 , using a coarse private grid followed by a local refinement around the best feasible candidates. For privacy bounds we use Theorem 3.1 together with Proposition 3 which gives Vn (r) = V̄n (r) with κ = √ √ exp(− 2/ℓ) and Lemma B.1 that provides the generic O( n)-bound for ∆n (r). We emphasise that this latter bound is not tailored to the present spatial model and is likely conservative: it controls the posterior mean sensitivity uniformly over all bounded-response datasets, without exploiting the smooth plume structure or the specific GP kernel form. Consequently, the privacy-constrained hyperparameter search may over-regularise the posterior, making the observed privacy-utility tradeoff more pronounced than what would be obtained from a sharper problem-specific sensitivity analysis. This is analogous to the 1D setting where tighter bounds on ∆n (r) lead to substantially less pessimistic privacy calibration. After hyperparameter selection, a single probability threshold C is chosen to maximise the mean IoU of the decision sets {x : pD (x) ≥ C} over the sampled population. We also evaluate the actual one-path private release:

Tomasz Maciazek

for each replicate, posterior sample paths are drawn on a grid, smoothed by a fixed Gaussian filter and thresholded to form released excursion sets. This smoothing is data-independent post-processing and therefore does not affect the privacy guarantee. We report the non-private posterior-probability excursion estimate {x : pD (x) ≥ C} and the IoU of the smoothed one-path private releases against the true excursion set. Table 4 and Figures 7–8 show that the privacy-constrained posterior can still yield useful released excursion sets, even when its posterior-probability calibration is substantially degraded. The effect of the privacy constraint is most visible in the BCE values: the private hyperparameters satisfying ε < 10 have much smaller effective dimension than the unconstrained posterior and this leads to a large relative increase in BCE, especially in the low-noise regime. However, the IoU degradation is considerably milder. At low noise the non-private posteriorprobability set is almost exact, so the private one-path releases necessarily show a visible loss relative to this very strong baseline. Nevertheless, the released boundaries remain concentrated around the true excursion boundary and are still highly informative. As the noise level increases, the unconstrained excursion set itself becomes less accurate and the additional cost of privacy becomes less pronounced in relative IoU terms. This transition is reflected both quantitatively in the table and qualitatively in the figures: for NSR = 0.11 the non-private boundary nearly coincides with the truth, whereas for NSR = 1 the non-private estimate is already visibly imperfect, making the smoothed private one-path releases a reasonable approximation despite the DP-induced randomisation. Mξ

0.1

0.3

0.5

NSR 0.11 0.43 1.00 ε, unconstrained 107 105 2· 103 deff , unconstrained 211.8 [210.8, 212.8] 96.63 [96.24, 97.00] 39.55 [39.42, 39.67] deff , ε < 10 5.64 [5.63, 5.64] 5.59 [5.58, 5.59] 5.47 [5.47, 5.48] relative BCE increase 527.5% [392.5%, 638.4%] 197.0% [153.6%, 240.5%] 91.5% [69.9%, 114.3%] IoU, unconstrained 0.979 [0.974, 0.983] 0.950 [0.939, 0.959] 0.908 [0.888, 0.927] 1-path IoU mean, ε < 10 0.824 [0.798, 0.850] 0.800 [0.772, 0.829] 0.776[0.743, 0.808] 1-path IoU SD, ε < 10 0.034 [0.029, 0.039] 0.039 [0.033, 0.045] 0.048 [0.041, 0.057] relative IoU gap 16.8% [13.8%, 19.4%] 15.8% [12.7%, 18.4%] 14.5% [10.4%, 18.8%] Table 4: The effect of varying the NSR set size on the released 2D excursion set utility at fixed training set size |D| = 400. Median [interquartile range] over 103 draws of (D, f∗ )-pairs across a range NSR values. The private candidate is the best average-BCE-selected setting satisfying ϵ < 10.

D.2

Greater London property sales study

London house-price excursion experiment. We also evaluate the method on a real spatial dataset derived from the HM Land Registry Price Paid Data. We restrict attention to Greater London leasehold-flat transactions in 2018 and join transaction postcodes to a coarse grid. Transactions are aggregated onto a hexagonal lattice of Greater London and each non-empty hexagon is assigned the median log sale price of the transactions falling in that cell. For a fixed threshold t, we define the centred response yhraw = mediani∈h log Pi − t, and the corresponding excursion label sh = 1{yhraw > 0}. The GP is fitted to the symmetrically clipped and rescaled response yh =

clip (yhraw , −B, B) ∈ [−1, 1], B

so that the response bound entering the privacy calculation is MY = 1. This transformation preserves the zero excursion threshold, apart from saturation of extreme values, since yhraw > 0 if and only if yh > 0. We use the same exponential-kernel posterior model as in the synthetic experiments, with hexagon-centre coordinates normalised to [0, 1]2 . For the non-private baseline hyperparameters (ℓ, r, σ) are selected by validation

Differential Privacy of Gaussian Process Posterior Sampling

Figure 7: Representative examples of the 2D excursion-set experiment at NSR = 1. Each row uses one independently sampled latent field f∗ and dataset D. Left: latent field f∗ and training data locations with the black curve indicating the true excursion boundary. Middle: non-private posterior-probability excursion estimate {x : pD (x) ≥ C} overlaid on the true excursion set shown in grey. Right: private one-path excursion releases obtained by drawing posterior sample paths, applying a smoothing post-processing step, and thresholding at the natural excursion level. The private (εL < 10, L = 1) one-path boundaries showcase the DP randomisation and associated IoU tradoff while still tracking the true excursion boundary well. Dataset size |D| = 400, δ = 0.001. binary cross-entropy (using separate transaction data from 2017), and the probability cutoff C defining the excursion estimate b C = {x : pD (x) ≥ C} Ω is selected by validation IoU. For the private version, the same selection procedure is constrained to hyperparameter settings satisfying the Rényi-DP-based privacy bound at the prescribed privacy level. The resulting private visualisations are produced by drawing sample paths from the DP-constrained posterior and thresholding them at zero. The resulting privacy guarantee should be interpreted at the level of the coarsened hexagon records used by the GP mechanism. Exact transaction locations are discarded before fitting: postcode-level coordinates are assigned to the coarse hexagon cells and the mechanism operates only on the resulting hexagon-level median responses. In Figure 4 the non-private probability boundary pD (x) = C is compared with three private one-path excursion boundaries. The private draws recover the dominant high-price region in central London, but tend to smooth out or miss smaller high-price islands outside central London, reflecting the additional randomisation and regularisation induced by the privacy constraint.

Tomasz Maciazek

Figure 8: A representative example of the 2D excursion-set experiment at NSR = 0.11. Left: latent field f∗ and training data locations with the black curve indicating the true excursion boundary. Middle: non-private posterior-probability excursion estimate {x : pD (x) ≥ C} overlaid on the true excursion set shown in grey. Right: private one-path excursion releases obtained by drawing posterior sample paths, applying a smoothing post-processing step, and thresholding at the natural excursion level. The private (εL < 10, L = 1) one-path boundaries showcase the DP randomisation and associated IoU tradoff while still tracking the true excursion boundary well. Dataset size |D| = 400, δ = 0.001.

Record · ID 282720 · SHA-256 f11f7ded1b7e7781
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.